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Notation 


Major  notation  is  listed  below.  Additional  characters  that  have  been  used  only  once  are  defined 
where  they  are  used. 

A  Flow  area 

Cp  Specific  heat 

Cd  Drag  coefficient 

Cvm  Virtual  mass  term  coefficient 

D  Diameter 

Dh  Hydraulic  diameter 

F  Force  term 

f  Friction  factor 

G  Mass  flux 

g  Gravitational  constant 

h  Specific  enthalpy 

hfg  Latent  heat  of  vaporization 

h  Convection  heat  transfer  coefficient 

K  Pressure  loss  coefficient 

k  Thermal  conductivity;  diamensionless  constant 

M  Molecular  mass  number 

m  Mass  fraction 

m  Mass  flow  rate 

P  Pressure 

Pi  Interfacial  perimeter 

q '  Thermal  load,  per  unit  length 
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q"  Heat  flux 

R  Universal  gas  constant 

Rcu  Radius  of  curvature 

Re  Reynolds  number 

S  Slip  ratio 

s  Specific  entropy 

T  Temperature 

t  Time 

U  Velocity 

V  Volume 

Wecr  Critical  Weber  number 

W  Power 

Xlm  Lockhart-Martirelli  parameter 
x  Quality;  mass  fraction 

Y  Mole  fraction 

z  Axial  coordinate 

Greek  Letters 

a  Void  fraction 

T  Evaporation  mass  rate,  per  unit  volume 

S  Passage  width 

e  Liquid  volume  fraction 

cr  g  Geometric  standard  deviation 


4 


e  Liquid  volume  fraction 

tjT  Turbine  isentropic  efficiency 

G  Angle  of  inclination  with  horizontal  plane 

p  Density 

Subscripts 

B  Bubble 

e  Equilibrium 

F  Frictional 

f  Saturated  liquis 

G  Vapor-noncondensable 

g  Saturated  vapor 

h  Homogeneous 

i  Interphase 

L  Liquid 

LS  Liquid  superficial 

n  Noncondensable 

s  Insentropic 

TP  Two-phase 

v  Vapor 

1,2  Inlet,  outlet 
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1. 


Introduction 


Fuel  cells  are  electrochemical  systems  that  directly  covert  chemical  energy  from  a  stream 
of  fossil  fuel  into  electricity,  by  combining  hydrogen  and  oxygen  with  the  aid  of  a  catalyst.  Fuel 
cells  do  not  apply  the  intermediate  step  of  converting  chemical  energy  into  heat  before  producing 
electricity,  and  thus  have  lower  irreversibility,  have  fewer  moving  parts,  and  do  not  generate 
many  of  the  harmful  combustion  products  typical  of  fossil  fuel  generators. 

Fuel  cells  typically  include  an  anode  and  a  cathode,  separated  by  an  electrolyte.  Several 
different  types  of  electrolytes  are  used,  including  the  proton  exchange  membrane,  concentrated 
phosphoric  acid,  molten  carbonate,  and  solid  ceramic  electrolyte. 

In  most  fuel  cells,  a  mixture  of  water  vapor  and  hydrogen  constitutes  the  anodic  fluid, 
while  the  cathodic  fluid  is  made  of  water  vapor  and  oxygen  (air).  Although  the  cathodic  fluid  is 
often  discarded,  recycling  of  water  content  of  the  cathodic  fluid  may  be  necessary  when  the  fuel 
cell  is  used  in  places  where  water  is  scarce.  For  recycling  of  the  water  content  of  the  cathdolic 
fluid,  water  should  first  be  condensed  in  order  to  separate  it  from  air.  The  condensate  must  then 
be  evaporated  before  flowing  into  the  anodic  flow  passages. 

The  design  of  the  evaporator  and  condenser  components  of  fuel  cell  systems,  however, 
requires  special  care,  due  to  their  potentially  unconventional  configurations  and  working 
conditions.  The  evaporator  may  include  small  (micro)  passages.  The  condenser,  on  the  other 
hand,  requires  to  be  specially  designed  due  to  the  very  high  concentrations  of  noncondensable 
(air)  in  the  cathodic  fluid. 

The  purpose  of  the  effort  reported  in  this  document  was  to  develop  engineering  analytical 
tools  for  the  design,  analysis,  interpretation,  and  extrapolation  of  experimental  data  to  be 
obtained  during  a  number  of  novel  evaporator  and  condenser  systems.  These  included  an 
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evaporator  using  a  heated  helicoidal  microchannel;  an  evaporator  using  a  large  number  of 
parallel  michrochannels  connected  at  their  two  ends  to  common  plenna;  an  evaporator  using  a 
large  number  of  parallel  microchannels  connected  at  their  two  ends  to  common  manifolds  and 
condensation  (due  to  spontaneous  fog  formation,  or  spray  assisted),  and  droplet  removal  in  a 
system  consisting  of  an  expander  and  a  cyclone  separator. 

In  what  follows,  the  basic  modeling  features  of  a  system  consisting  of  a  single  heated 
michrochannel  evaporator  are  first  briefly  discussed  in  Chapter  2.  The  models  dealing  with 
multi-channel  evaporators  are  then  described  in  Chapter  3.  In  Chapter  4  the  models  developed 
for  the  condenser  and  moisture  separation  systems  are  described.  Draft  users  manuals  for  several 
computer  codes  that  have  been  completed  in  this  project  are  then  provided  as  Appendices. 

Two  other  major  tasks  were  also  carried  out  in  this  project,  which  are  not  reported  in  this 
document,  since  they  were  not  fully  completed.  These  included  the  development  of  a  surface 
condenser  model  based  on  the  porous-media  representation  of  the  secondary  side;  and  a  generic 
one-dimensional  model  for  the  simultaneous  fog  formation/condensation  on  droplets  and 
droplets  removal  due  to  droplet  impaction  on  flow  passage  surfaces.  The  former  task  was  under 
taken  early  during  the  project,  but  was  halted  when  a  different  condenser  design  was  selected. 
The  latter  task  was  undertaken  with  the  intention  of  developing  a  simple  but  flexible  and  general- 
purpose  model  for  the  simulation  of  condensation  and  droplet  removal  in  flow  passages  with 
arbitrary  geometry. 
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2.0  Single-Channel  Evaporator  Models 
2.1  General  Remarks 

Two-phase  flow  and  boiling  in  microchanels  and  small  passages  is  of  great  interest  and 
has  been  extensively  studied  recently  (Ghiaasiaan  and  Abdel-Khalik,  2000).  Micropassages  and 
networks  of  microchannels  can  provide  extremely  large  surface-to-volume  ratios,  and  can  thus 
dispose  of  extremely  large  volumetric  thermal  loads.  Most  of  the  applications  involving  cooling 
systems  are  based  on  subcooled  forced  convection  and  nucleate  boiling,  where  large  heat  fluxes 
can  be  sustained  with  relatively  low  wall  temperatures.  For  such  applications  the  onset  of 
nucleate  boiling  and  onset  of  flow  instability,  and  critical  heat  flux  are  important  operational  and 
safety  thresholds.  Other  issues  of  interest  include  pressure  drop,  two-phase  flow  patterns  and 
their  characteristics,  and  void  fraction.  Only  a  brief  review  of  issues  relevant  to  narrow  channels, 
which  are  of  interest  to  this  project,  will  be  provided. 

Investigations  recently  published,  that  deal  with  two-phase  flow  in  micropassages  include 
Ali  and  Kawaji  (1991),  Ale  et  al.  (1993),  Mishima  et  al.  (1993),  Wilmarth  and  Ishii  (1994),  and 
Fourar  and  Bories  (1995).  Depending  on  the  liquid  and  gas  (vapor)  flow  rates,  several  flow 
patterns  may  develop.  For  vertical  channels,  the  void  fraction  can  be  obtained  using  the  drift 
flux  model  (DFM),  and  Mishima  et  al.  (1993)  has  proposed  model  constants  applicable  to  thins 
passages. 

An  evidently  important  issue  related  to  two-phase  flow  and  boiling  in  thin  passages  is 
pressure  drop.  The  aforementioned  authors,  with  the  exception  of  Wilmarth  and  Ishii  (1994),  all 
studied  pressure  drop  as  well.  For  laminar  and  turbulent  single-phase  flow,  Jones  (1974)  has 
developed  a  simple  and  accurate  method  for  the  prediction  of  friction  factors.  A  useful 
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correlation  applicable  to  very  narrow  channels  with  rough  surfaces  has  also  been  derived  by  John 
et  al.  (1988).  Most  of  the  authors  who  have  measured  two-phase  pressure  drop,  have  attempted 
to  modify  the  constants  in  the  correlation  of  Chisholm  and  Laird  (1958).  The  resulting 
modifications,  however,  are  often  complicated  and  appear  to  be  system-specific  (Mishima  et  al., 
(1993).  Recently,  Ekberg  (1999)  measured  pressure  drops  in  horizontal  thin  annuli.  A  simple 
pressure  drop  correlation  previously  derived  by  Beattie  and  Whalley  (1982)  was  subsequently 
found  to  well-predict  the  latter  data.  The  correlation  of  Beattie  and  Whalley  may  thus  be 
considered  as  an  accurate  and  simple  method  for  two-phase  pressure  drop  calculation  in  thin 
micropassages. 

Flow  and  heat  transfer  in  helicoidal  passages  have  been  studied  rather  extensively  in  the 
past,  since  helicoidal  channels  are  used  in  compact  heat  exchangers,  chemical  reactors, 
refrigeration  and  cryogemic  systems.  Recent  studies  of  single-phase  flow  include  Lin  and 
Ebadian  (1997)  and  Li,  Lin  and  Ebadian  (1998).  Friction  factor  and  heat  transfer  coefficient  are 
both  enhanced  in  comparison  with  straight  channels.  Most  of  the  published  studies  deal  with 
circular  channels,  however,  and  non-circular  channels  have  received  little  attention  (Shah  and 
Joshi,  1987). 

For  two-phase  flow,  the  coil  orientation  and  helix  angle  are  very  important  parameters. 
For  vertically  oriented  helicoidal  channels  with  small  angle  the  flow  regimes  appear  to  have 
much  in  common  with  those  observed  in  horizontal,  straight  channels  (Kaji  et  al.,  1984). 
Published  studies  dealing  with  two-phase  flow  regimes  and  pressure  drop  include  Baneijee  et  al., 
(1969),  Saxena  et  al.  (1990),  Awwad  et  al.  (1994  a,b;  1995),  Xin  et  al.  (1990),  Kang  et  al. 
(2000).  Unfortunately,  however,  the  past  investigations  were  primarily  concerned  with  circular 
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channels,  and  data  and  analysis  dealing  with  two-phase  flow  in  non-circular  helicoidal  flow 
passages  are  scarce. 

2.2  Conservation  Equations 

One-dimensional  conservation  equations  representing  flow  in  a  heated  channel  with 
uniform  cross-section  are  presented  in  this  section. 

A  schematic  of  a  heated  channel  is  displayed  in  Fig.  1.  In  general,  the  flow  field  can  be 
divided  into  at  least  four  zones;  single-phase  subcooled  liquid  near  inlet  (Zone  1);  a  single-phase 
superheated  vapor  near  the  exit  (Zone  4);  and  boiling  two-phase  flow  in-between.  The  boiling 
two-phase  flow  region  is  also  divided  into  two  different  zones  here  for  convenience.  In  Zone  3 
all  boiling  regimes  are  considered.  High  quality  (high  void  fraction)  flow,  representing  the 
evaporation  of  an  extremely  thin  liquid  film  in  annular  flow  regime,  however,  is  dealt  with 
separately  in  Zone  4,  since,  for  fast  and  efficient  numerical  solution,  the  two-phase  flow 
conservation  equations  need  to  be  presented  slightly  differently  in  the  latter  flow  regime. 

The  following  general  assumptions  are  made: 

1 .  The  flow  is  one  dimensional  and  steady-state 

2.  Channel  cross-section  is  uniform 

3.  Where  vapor  and  liquid  are  both  present,  thermal  equilibrium  between  the  two 
phases  exits 

4.  Wall  heat  flux  is  known 

5.  The  fluid  is  single-component 

2.2.1  Conservation  Equations  for  Single-Phase  Flow  (Zone  1  and  4) 
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The  liquid  -phase  mass  and  momentum  conservation  equations  are: 

(!-«)]■ =-r  0) 
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where  T  is  the  rate  of  mass  evaporation  per  unit  length  of  the  channel.  Elimination  of  T 
between  the  above  two  equations  gives: 
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Equations  (4),  (5),  (6),  and  (9)  are  four  coupled  ordinary  differential  equation  with  P ,  Uf,  Ug 
and  a  as  their  state  variables. 

2.2.3  Conservation  Equation  for  the  Film  Evaporation  Regime  (Zone  3) 

Although  the  conservation  equations  of  previous  section  are  directly  appreciable  here,  the 
numerical  solution  of  the  equations  is  somewhat  problematic  in  this  regime.  The  difficulty  is  due 

dec 

to  the  fact  that  in  this  regime  a  « 1  and  —  «  o ,  and  these  render  the  set  of  ordinary  differential 

dz 

equations  (ODEs)  stiff.  Furthermore,  since  the  two-phase  equations  must  in  principle  be 
integrated  until  x  =  \  and  a  =  1.0  (i.e.,  complete  evaporation  of  the  annular  film),  the  set  of 
ODEs  becomes  stiffer  as  the  integration  proceeds.  The  liquid-phase  momentum  equation  (Eqn. 
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9)  is  the  main  source  of  difficulty  since  it  represents  the  acceleration  of  a  fluid  with  essentially 
zero  inertia. 

To  avoid  this  difficulty,  the  two-fluid  scheme  of  the  previous  section  is  simplified  by 
dropping  Eqn.  (9),  (i.e.  the  liquid  phasic  momentum  equation).  Instead,  the  following  slip  ratio 
expression  of  Zivi  (1968)  is  used: 


S  =  Ug/Uf=(pf/ps)m  (9-a) 

The  above  expression  has  been  derived  based  on  the  assumption  of  minimum  entropy 
generation  in  annular  flow. 

Using  Eqn.  (9a),  Ug  is  eliminated  from  Eqn.  (4-6),  leading  to  three  ODEs,  with  P,  U f 
and  a  as  the  state  variables.  The  void  fraction,  furthermore,  is  eliminated  in  favor  of 
(equilibrium)  quality,  .  The  resulting  equations,  which  are  given  below,  can  be  integrated 
without  difficulty. 


1  <K 

xe  dz 


M 
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dP  1  dU 
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dz  Ug  dz 
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It  should  be  noted  that,  deriving  Eqn.  (10-a),  —  «  o  has  bee  assumed. 
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2.3  Flow  Regime  Maps  and  Interfacial  Force  Terms 
2.3.1  Straight  Channels 

Two-phase  flow  regimes  in  narrow,  rectangular  channels  are  complicated,  since  velocity 
slip  readily  occurs  in  such  channels.  A  good  review  can  be  found  in  Ghiaasiaan  and  Abdel- 
Khalik  (2001). 

A  good  representative  of  previously  published  studies  is  that  of  Ali  and  Kawaji  (1991). 
The  observed  regimes  in  their  channels  generally  include  bubbly  (small  bubbles);  intermittent, 
and  annular.  The  intermittent  regime  includes  slug,  and  chum  flow  patterns.  Rivulet  flow  can 
also  occur  at  very  low  liquid  superficial  velocities.  Complicated  flow  patterns  were  also  reported 
by  Fourar  and  Bories  (1995).  Published  relevant  data  have  significant  disagreement,  and  no 
reliable  flow  regime  map  is  currently  available  (Ghiaasiaan  and  Abdel-Khalik,  2001).  Therefore, 
in  this  study,  a  simple  but  well-tested  flow  regime  transition  model  will  be  used.  This  simple 
method  has  been  successfully  used  in  the  past  for  two-phase  flow  in  channels  undergoing  rapid 
evaporation.  (Ghiaasiaan  and  Geng,  1997).  Accordingly,  we  assume: 

Bubbly  for  a  <  0.2  (11) 

Annular  for  a  >  0.80  (12) 

Intermittent  for  0.2  <  a  <  0.75  (13) 

In  bubbly  regime,  we  assume  bubbles  with  diameters: 


Db  =min 


MUe-Uj 


We  — 

2  "  Kcr  ’  r, 


(14) 


where  S  =  width  of  the  channel,  and 


Equation  (14)  assumes  that  aerodynamic  dispersion  is  responsible  for  bubble  breakup.  The  void 
fraction  and  number  density  of  bubbles  in  bubbly  flow  regime  are  related  according  to 


TV'  =  — — —  (16) 

— D  3 


Furthermore,  the  interfacial  are  a  concentration  can  be  found  from: 

^  =  NttDb2  (17) 

A 

The  interfacial  force  per  unit  mixture  volume  can  now  be  found  from: 

F‘=\a%p'\u‘-u'  I  (18) 

where: 

Q,  =^-(1  +  0.15  Re/687)  (19) 

Re5  v 

RcB=pf\Ug-Uf\DB/Mf  (20) 

In  the  annular  flow  regime,  it  can  be  shown  that 

F,^\Us-Uf\  (21) 

where  it  is  assumed  that 

Cfl  =  0.05  (22) 

In  the  intermittent  regime  (see  Eqn.  (13)),  Fi  is  found  by  interpolation  between  the 
bubbly  and  annular  regime  values  at  a  =  0.25  and  a  =  0.8 ,  respectively,  using  the  void  fraction 
as  the  interpolation  parameter. 
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The  virtual  mass  force  term,  the  last  term  on  the  right  side  of  Eqns.  (8)  or  (9),  involves 
the  constant  Cm .  This  coefficient  is  represented  as: 

cm  =  a(l-a)  [pf{\-a)+pgct]  (23) 

2.3.2  Helicoidal  Channels 

The  above  models,  strictly  speaking,  are  appreciable  to  straight  channels  and  may  need  to 
be  modified  before  they  can  be  applied  to  helicoidal  channels.  Little  is  known  about  the  flow 
regimes  in  helicoidal  channels  that  have  rectangular  cross-sections,  however.  A  brief  discussion 
follows. 

For  vertically-oriented  helicoidals,  when  the  helix  angle  is  very  small,  and  provided  that 
Dh  /  Rcu  «1  (where  Reu  represents  the  coil  radius  of  curvature),  the  flow  regimes  can  be 
expected  to  be  similar  to  those  occurring  in  straight  channels,  and  commonly-applied  methods 
for  two-phase  flow  pressure  drop  can  also  be  used  (Baneijee  et  al.,  1969).  Kaji  et  al.  (1984)  also 
noted  similarity  between  straight,  horizontal  channels  and  vertical  helicoidal  channels  with  small 
helix  angles,  with  respect  to  flow  regimes.  A  good  experimental  study  aimed  at  the 
measurement  and  correlation  of  two-phase  friction  factors  in  vertical  helicoidal  channels  with 
circular  cross-section  with  helix  angles  in  the  0.5°  □  10°,  has  been  reported  by  Xin  et  al.  (1996). 
The  developed  empirical  correlation  for  friction  factor  in  the  latter  study  of  course  may  not  be 
directly  applicable  to  noncircular  channels. 

Awwad  et  al.  (1995)  measured  and  correlated  the  air-water  two-phase  pressure  drop  in  a 
horizontally-oriented  helicoidal  tube,  for  helix  angles  in  the  1.0-  20°  range,  and  pipe  diameters  in 
the  330  -  670mm  range.  The  single-phase  friction  factors  in  their  experiments  agreed  well  with 
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the  correlations  of  Manlapaz  and  Churchill  (1980)  for  laminar  flow,  and  Ito  (1959)  for  turbulent 
flow.  They  correlated  their  two-phase  flow  pressure  drop  data  as: 


<t>L  = 


\+Xlm 


cf;a 


r  12  1  V/2 

1  + - +- 


v 


^LM  XLM  J 


(24) 


where  XLM  is  the  Martinelli  factor,  and  n  =  0.576 ,  and  C  =  7.79  for  Fd  <  0.3 ;  while  «  =  1.3  and 
C  =  13.56  otherwise.  The  dimensionless  number  Fd  is  defined  as: 

Fd=Fk(DnRjA  (25) 

gD 

where  ULS  is  the  liquid-phase  superficial  velocity  and  D  is  the  (circular)  channel  diameter. 

For  horizontal  helicoidal  channels,  due  to  the  apparent  lack  of  better  and  more  relevant 
data  and  correlations,  Eqn.  (24)  can  be  used  for  rectangular  channels,  with  D  replaced  with  Du , 
albeit  that  the  correlation  has  been  developed  originally  for  circular  channels.  Should 
horizontally-oriented  helicoidal  channels  be  of  interest,  furthermore,  the  flow  regime  model  and 
the  interfacial  force  terms  described  in  the  previous  section  will  be  unlikely  to  be  accurate,  and 
may  need  modification. 

For  the  vertical  helicoidal  configuration,  the  methodology  described  in  the  previous 
section  is  applicable  when  the  helix  angle  is  very  small,  and  for  circular  cross-section  channels. 
Major  modifications  may  be  needed  otherwise,  however. 


2.4  Wall  Friction 

Since  rectangular  channels  are  of  interest,  only  models  relevant  to  this  geometry  are  used 
at  this  time.  For  single-phase  flow,  and  assuming  a  straight  channel,  laminar-to-turbulent  flow 
regime  transition  is  assumed  to  occur  when  Re>2,300,  where  the  Reynolds  number  is  defined 
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based  on  the  hydraulic  diameter.  For  laminar  flow,  the  correlation  of  Shah  and  London  (1987) 
(See  Kakac  et  al.,  1987)  is  used,  whereby  the  D’Arcy  friction  factor  /  is  obtained  from: 

/Re  =  (24)|\-1.3553a*+L9467a'2  -1.7012a*3  +  0.9564a*4  -0.2537a*5  (26) 

where  a*  is  the  ratio  between  smaller  and  larger  sides  of  the  rectangular  cross-section.  For 
turbulent  flow,  the  correlation  of  Shah  and  Bhatti  (1987),  (KaKac  et  al.,  1987)  is  used,  where  by: 

F  =  (1.0875-0.1125  a)fc  (27) 

where  fc  is  the  circular-channel  correlation  of  Techo  et  al  (KaKac  et  al.,  1987): 
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The  two-phase  Reynolds  number  defined  in  Eqn.(28)  is  now  used  in  single-phase  friction  factor 
correlations  in  order  to  derive  the  two-phase  flow  friction  factor,  fTP .  The  frictional  pressure 
gradient  in  the  channel  is  than  obtained  from: 

(-dP/dz)F=f„-f-  (34) 

^ Ph 

It  should  be  mentioned  that  the  correlation  of  Beatttie  and  Whalley  is  based  on  circular  channel 
data,  Mid  uses  the  well-known  Colebrook- White  correlation  for  turbulent  single-phase  friction 
factor. 

It  should  also  be  emphasized  that  the  above  methods  may  need  to  be  modified  for 
application  to  helicoidal  channels. 

2.5  Propetties 

Thermodynamic  and  transport  properties  of  subcooled  and  saturated  water,  and  saturated 
and  superheated  steam,  are  calculated  using  the  routine  WSPROP,  previously  developed  and 
used  by  author  and  co-workers  (Ghiaasiaan  et  al.,  1995).  WSPROP  contains  several  explicit 
approximating  expressions  for  calculating  liquid  water  and  steam  thermodynamic  properties 
borrowed  from  Kuck  (1982).  The  partial  derivatives  of  water  and  steam  thermodynamic 
properties  are  also  calculated  in  WSPROP  by  differentiating  the  relevant  explicit  approximating 
expressions  (Kuck,  1982).  Water  and  steam  thermal  conductivity  and  viscosity  are  obtained  by 
interpolation  on  relevant  property  tables  of  Grigull  et  al.,  (1984).  The  interpolation  algorithms 
are  based  on  Press  et  al.  (1988). 
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2.6  Numerical  Solution  Method 


The  conservation  equations  were  expanded  using  appropriate  thermodynamic  relations 
and  cast  in  the  form  of  coupled  sets  of  ordinary  differential  equations  (ODEs),  using  the 
methodology  described  elsewhere  (Ghiaasiaan  and  Geng,  1997;  Ghiaasiaan  et  al.,  1995).  The 
details  can  be  found  in  the  latter  references  and  will  not  repeated  here.  The  resulting  sets  of 
ODES  have  the  following  generic  form: 
dY 

A——  =  C  (35) 

dz 

when  Y  is  a  column  vector  with  N  elements  the  state  variables,  with  //representing  the  number 
of  equation,  and  A  is  an  NxN  coefficient  matrix.  The  column  vector  C  is  an  N-element  vector. 

For  subcooled  single-phase  liquid  or  (superheated)  vapor  flow,  given  the  fact  that  G  = 
const,  (see  Eqns.  (1-3), 

Y  =  (hk,Pf  (36) 

with  k  =  L  and  v  for  liquid  and  vapor  respectively.  For  the  two-phase  flow  region. 

Y  =  (P,Ug,Upa)T  (37) 

Finally,  for  the  film  evaporation  regime: 

Y  =  {P,Ug,xe)T  (38) 

where  xe  is  the  local  equilibrium  quality. 

The  above  sets  of  ODEs  are  numerically  integrated,  in  tandem,  using  the  DVODE 
equation  solver  package  which  is  similar  to  the  equation  solves  described  by  Hindmarsh  (1980). 
This  numerical  package  uses  a  robust  and  efficient  variable-step  and  variable-order  numerical 
method  that  is  particularly  useful  for  the  solution  of  stiff  equation  systems. 
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3.0  Multi-Channel  Evaporator  Models 

3.1  General  Remarks 

Multi-channel  evaporators,  consisting  of  a  large  number  of  parallel  and  geometrically- 
identical  straight  heated  channels  connected  to  two  common  plena  or  manifolds  at  their  two 
ends,  have  been  considered  as  possible  designs  for  the  fuel  cells  of  interest.  The  models  briefly 
described  below  have  been  developed  for  such  multi-channel  evaporator  systems. 

It  should  be  noted  that  the  multi-channel  models  described  here  utilize  the  single-channel 
models  described  in  chapter  2  of  this  report  for  analyzing  the  flow  and  heat  transfer  phenomena 
in  each  individual  channel.  The  crucial  calculations  in  the  multi-channel  model  deal  with  the 
correct  implementation  of  the  inlet  and  exit  channel  boundary  conditions,  and  this  can  only  be 
done  by  carefully-designed  iterative  schemes. 

3.2  Multi-Channel  Evaporator  with  Inlet  and  Exit  Plena 

The  modeling  of  this  type  of  multi-channel  evaporators  is  actually  relatively  simple, 
although  the  iterative  solution  procedures  needed  for  this  purpose  are  computationally  expensive. 
The  simplicity  results  from  the  fact  that  the  parallel  channels  all  have  common  inlet  and  exit 
pressures. 

The  following  assumptions  are  made: 

1 .  The  parallel  channels  have  the  same  length  and  thermal  loads. 

2.  The  channels  can  be  divided  into  a  finite  number  of  groups.  In  each 
group,  the  channels  are  geometrically  identical. 

3.  Inlet  pressure  and  fluid  conditions  are  the  same  for  all  channels. 
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4.  The  channels  all  open  to  the  same  plenum,  and  therefore  have  the  same 
exit  pressure. 

Accordingly,  an  algorithm  was  developed  that  can  be  used  to  distribute  a  known,  total 
mass  flow  rate  among  a  given  set  of  parallel  channels.  The  algorithm,  briefly,  is  as  follows. 

The  average  mass  flux  in  the  channels  is  first  calculated.  Using  this  mass  flux  as  the  first 
guess,  the  conservation  equations  described  in  Chapter  2  are  solved  for  a  channel  from  each 
channel  group,  leading  to  the  calculation  of  channel  exit  pressures.  The  pressure  drops  of  all 
channel  groups  are  thus  obtained  and  are  averaged.  The  pressure  drop  in  each  group  is  then 
compared  with  the  aforementioned  average  pressure  drop.  The  mass  flux  for  each  tube  group  is 
then  re-adjusted  by  multiplying  it  by  the  ratio  of  channel  exit  pressure  calculated  in  the  previous 
iteration  and  the  average  exit  pressure  of  all  channels.  The  entire  procedure  is  then  repeated. 
Convergence  is  obtained  when  the  differences  among  the  calculated  exit  channel  pressures  are 
all  within  a  specified  tolerance  value. 

The  computer  program  that  performs  this  type  of  calculations  (EYAP10D)  has  been 
successfully  tested  for  an  imaginary  system  consisting  of  three  groups  of  tubes,  each  group 
consisting  of  five  geometrically-identical  channels. 

However,  the  algorithm  at  its  present  form  converges  relatively  slowly,  and  the 
aforementioned  computer  program  is  consequently  computationally  expensive. 

3.3  Multi-Channel  Evaporator  with  Inlet  and  Exit  Manifolds 
3.3.1  Mathematical  Model 

The  system  to  be  modeled  is  depicted  in  Figure  2.  The  geometry  clearly  shows  that  the 
inlet  and  exit  conditions  for  the  individual  channels  will  be  different,  and  the  distribution  of  the 
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fluid  mass  among  the  channels  will  be  non-uniform.  Once  again,  a  method  must  be  devised  to 
iteratively  distribute  the  mass  among  the  channels  such  that  all  boundary  and  interface  conditions 
are  correctly  satisfied. 

The  following  assumptions  are  made. 

1.  The  flow  in  all  channels,  and  along  both  manifolds,  is  one-dimensional  and 


steady-state. 

2.  The  fluid  at  the  exit  of  all  the  channels  is  pure  saturated  or  superheated  vapor. 

3.  At  inlet  to  the  inlet  plenum,  the  fluid  is  subcooled  liquid. 

4.  For  pressure  drop  calculation  purposes,  the  vapor  throughout  the  exit  plenum  can 
be  assumed  to  be  pure  saturated  vapor. 

5.  The  channels  are  geometrically  identical. 

6.  The  thermal  load  of  all  the  channels  are  the  same. 

In  what  follows,  the  fluid  momentum  equations  are  first  derived  using  the  vector 
momentum  method.  The  algorithm  used  for  the  iterative  solution  of  the  equation  is  then  briefly 
described. 

We  will  start  with  the  control  volume  depicted  in  Fig.  3,  which  represents  the  segment  of 
the  inlet  manifold  between  the  manifold’s  inlet,  and  the  inlet  to  the  first  heated  channel.  Assume 

that  the  static  pressure  at  inlet  to  the  inlet  plenum  is  P*  .  This  pressure  evidently  depends  on  the 
details  of  the  flow  field  upstream  from  the  inlet  plenum,  and  is  assumed  known. 

Applying  the  momentum  conservation  principle,  in  the  direction  of  the  inlet  manifold,  to 
the  depicted  control  volume,  one  gets: 
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The  parameters  in  the  above  equations  have  the  following  definitions: 

Amini  =  Inlet  manifold  cross-sectional  area 

Al  =  Channel  #1  cross-sectional  area 

Dh  m  in  =  Inlet  manifold  hydraulic  diameter 

Kent  =  Manifold  entrance  pressure  loss  coefficient 

KTent  =  Pressure  loss  coefficient  due  to  entrance  in  a  T-junction 

f.  =  D’Arcy  friction  factor  for  the  inlet  manifold 

Other  parameters,  including  pressures  and  mass  flow  rates  are  defined  in  Fig.  3. 

Now,  consider  the  control  volume  depicted  in  Fig.  4,  which  represents  the  segment  of  the  inlet 
manifold  between  channels  j- 1  and  j .  For  this  control  volume  one  can  write 
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The  definitions  of  parameters  in  the  above  equation  are: 


thj  =  Mass  flow  rate  in  the  jth  channel 


Area  of  the  y'th  channel 


fj  = 


Friction  factor  in  the  depicted  segment  of  the  manifold 


S  =  The  spacing  between  two  adjusted  channels 

Flow  rates  and  pressures  are  defined  in  Fig.  4.  Figure  4,  and  Eqns.(40-42)  in  fact  apply 
to  channels  2  through  N ,  with  N  representing  the  total  number  of  channels. 

Similar  equation  sets  can  be  set  up  for  the  exit  manifold.  Figure  5  depicts  the  control 
volume  at  the  exit  of  the  exit  manifold.  For  this  control  volume  one  can  write: 
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where: 


Kexit  =  Manifold  exit  loss  coefficient 


Kt  .  -  Loss  coefficient  for  exit  T junction 
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Area  of  the  exit  manifold 


An,o  = 

Dhmo  =  Hydraulic  diameter  of  the  exit  manifold 
Likewise,  for  the  control  volume  depicted  in  Fig.  6,  we  can  write: 
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3.3.2  Solution  Method 

The  equations  presented  in  the  previous  section,  when  written  for  all  the  control  volumes 
in  the  inlet  and  exit  manifolds,  form  a  set  of  algebraic  equations  which  would  be  closed  for  the 
conditions  to  be  described  shortly,  where  the  single-channel  models  described  in  Chapter  2  are 
included.  In  other  words,  with  geometric  characteristics  of  a  particular  channel  as  well  as  its 
thermal  load  and  inlet  fluid  enthalpy  all  fixed,  the  entire  single-channel  model  can  be  thought  of 
as  a  highly  non-linear  algebraic  equation  in  the  generic  form: 

f(PinJ>PouU>™j)  =  0  (5°) 

The  above  equations,  including  those  represented  by  Eqn.(49)  for  j  =  1,2,  —N,  constitute  a 
closed  set  of  equations,  with  unknown  that  include:  m} ,  PmJnJPmt0j ,  PinJ ,  Pou,j ,  etc..  Furthermore, 

two  of  the  following  three  quantities  must  be  known:  P*n,P0Ut,and  mlot . 
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An  iterative  solution  of  there  equation  is  evidently  difficult  and  time  consuming.  After 
trying  a  number  of  different  methods,  the  following  novel  method  was  eventually  found  to  be  the 
most  efficient. 

We  consider  the  conditions  where  P'n  and  mwt  are  assumed  known.  The  following 
procedure  is  then  used.  (Note  that ,  if  P*„  or  mtot  is  the  unknown  the  method  is  still  valid,  only 
more  iteration,  for  example  on  P*n  when  it  is  an  unknown  are  needed). 

1.  Knowing  the  total  mass  flow  rate  and  inlet/exit  pressures  for  the  manifolds,  the  single¬ 

channel 

model  is  applied  to  generate  a  table  with  about  60  entries,  providing  predicted  values  of 
channel  exit  pressure  given  inlet  channel  pressure  and  channel  mass  flow  rate.  These 
runs  are  not  iterative  and  can  be  performed  with  relatively  little  computation.  The  entries 
in  this  table,  when  used  along  with  an  interpolation  routine,  in  fact  represent  the 
numerical  solution  for  Eqn.  (49). 

2.  The  conservation  equation  for  the  first  channel,  and  the  segments  of  the  inlet  and  exit 
manifolds  that  are  located  between  the  first  channel  and  their  entrance  and  exit, 
respectively,  are  iteratively  solved.  These  equations  includes  Eqn.  (38-40),  (44,45),  and 
(49)  with  y=l.  (Note  that  the  latter  equation  is  replaced  by  an  interpolation  routine  that 
interpolates  in  the  afore-mentioned  table).  The  equation  set  is  in  fact  closed,  and  would 

lead  to  the  calculation  of  mx ,  P*ul ,  mm  in  l ,  PmM ,  rhmoX ,  Pm  o>1 . 

3.  The  solution  in  step  2  above  provides  the  boundary  conditions  for  the  system  consisting 
of  the  channel  #2,  and  the  segments  of  the  inlet  and  exit  manifolds  that  are  located 
between  the  first  and  second  channels.  The  relevant  equations  included  Eqns.  (41-43), 
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and  (46-48),  and  (49),  and  their  iterative  solution  leads  to  the  calculation  of  m2 , 


P,nM,X  ’  <JnJ  ’  Pm,oJ  >  mm,oJ  >  Pi„J -  Pou<J  ^d  tttj  fOI  }=2. 


4.  Step  3  is  repeated  for  the  third  channel,  and  so  on. 

The  solution  would  end  at  the  last  channel,  for  which  the  predicted  flow  rate 
should  match  what  is  left  from  an  overall  mass  balance  for  the  entire  system.  The 
aforementioned  algorithm  has  been  developed,  and  tested.  Numerical  tests  show  that  the 
algorithm  converges  well.  Preliminary  application  of  the  model  to  a  test  data  previously 
provided  to  us  by  Tyndall,  however,  showed  that  for  the  given  boundary  conditions  the  total 
mass  flow  rate  could  have  been  carried  through  by  only  about  half  the  channels.  This 
indicates  that  we  may  be  under  estimating  the  pressure  drops  in  the  system.  Underestimation 
of  singular  pressure  losses  associated  with  the  numerous  bends  and  turns  in  the  system  is  of 
course  a  potential  main  contributor.  More  accurate  representation  of  the  entire  system 
geometry  is  evidently  needed. 
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Fig.4  The  Control  Volume  Between 
Entrances  to  Channels  j-1  and  j 
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Fig.5  The  Control  Volume  at  the  Exit  of 
the  Exit  Manifold 
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4.0  Condenser  and  Moisture  Separation  Systems 

4.1  General  Remarks 

Efficient  and  innovative  condenser  designs  are  needed  for  the  fuel  cells  of  interest  in  this 
study,  since  the  vapor  pressure  in  the  vapor  noncondensable  mixture  that  should  be  used 
for  condensation  is  typically  low.  Surface  condensers,  as  is  well-known,  are  inefficient 
when  large  concentrations  of  noncondensables  are  present.  For  this  reason,  the  condenser 
systems  of  interest  in  the  study  consist  of : 

1) .  An  expander,  in  which  the  vapor-noncondensable  is  depressurized.  Partial 

condensation  of  the  vapor  phase  then  results,  either  as  a  result  of  fog  formation,  or 
as  a  result  of  condensation  on  micro  droplets  of  water  that  have  intentionally  been 
sprayed  into  the  flown  field. 

2) .  A  cyclone-like  separator,  in  which  the  droplet-laden  gas-vapor  mixture  flows 

through  a  tortuous  passage.  The  highly  curved  streamlines  that  are  imposed  on  the 
mixture  result  in  the  inertial  impaction  of  the  entrained  droplets  on  the  walls  of  the 
flow  passage. 

In  what  follows,  engineering  models  are  briefly  described,  that  have  been  developed  for  the 
simulation  of  a  system  consisting  of  an  expander  with  known  inlet  and  exit  pressures,  and  a 
cyclone-type  separator,  that  are  configured  in  series. 

4.2  Expander  Models 

Four  types  of  expanders  are  modeled.  Expander  Types  1  and  4  consist  of  a  system,  with 
or  without  moving  parts,  with  known  inlet  and  outlet  pressures.  Expander  Types  2  and  3,  on  the 
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other  hand,  are  turbine  expanders  with  known  input  pressures,  known  power  outputs,  and  known 
isentropic  efficiencies.  For  these  expanders  the  exit  pressure  is  actually  calculated  using 
thermodynamic  relations.  In  Expander  Types  1  and  3,  furthermore,  droplets  with  a  log-normal 
size  distribution  at  a  constant  rate  are  assumed  to  be  injected  into  the  vapor-noncondensable 
mixture  at  inlet,  in  order  to  act  as  condensation  sites.  In  Expander  Types  2  and  4,  however, 
spontaneous  fog  formation  is  assumed  to  take  place. 

Common  to  all  of  the  afore-mentioned  expander  types  are  the  following  assumptions: 

1 .  The  system  is  at  steady-state. 

2.  The  vapor-noncondensable  mixture  and  its  entrained  liquid  droplets  are 
everywhere  at  thermodynamic  equilibrium. 

3.  The  gas  and  liquid  droplets  everywhere  have  the  same  velocity. 

4.  The  liquid  condensate  is  impermeable  to  the  noncondensable. 

5.  The  noncondensable  is  an  ideal  gas 

6.  Except  for  the  inlet  and  outlet,  no  mass  transfer  through  the  walls  of  the  system 
takes  places 

7.  Droplet  break-up  and  coalescence  are  neglected. 

4.2.1  Expander  Types  1  and  4 

First  consider  the  Expander  Types  1  or  4.  The  schematic  in  Figure  7  represents  the 
modeled  system.  For  this  system  the  following  balance  equation  can  be  written: 

•  Total  Mass  balances: 

AUX  [(!-£,)  A?,2+Gl  P/>l]  =  ^2t/2[(1-e2)/7G.2+G2  P/.2]  (51) 
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where  A  represents  the  flow  area,  U  represents  velocity,  and  e  represents  the  in-situ  volume 
fraction  occupied  by  liquid  droplets.  Subscripts  1  and  2  represent  the  inlet  and  exit,  respectively. 
Also,  subscript  /  represents  the  liquid  phase  while  subscript  G  represents  the  vapor 
noncondensable  mixture. 

•  Noncondensable  Mass  Balance: 


[0~~  €  )  Pc,\Xn,\  ]  =  Ar^J 2  [(l—  e2  )  Pc,2Xn,2  ]  (**2) 

where  xn  in  this  and  the  forthcoming  equations  represents  the  mass  fraction  of  the 
noncondensable  in  the  vapor-noncondensable  mixture. 

•  Mixture  Energy  Balance: 


f  l 

(l—  ei )  Pg,\  ^g,i  ■*"  T  ^ i 
V  2  J 


AXUX 
+Q-W  =  A2U2 


\  ( 
+  Pf,\ 


Kt+\u' 


(l-e2) 


V 


KiPui' 

**  J 


+  e2  P/,2  ^/,2  +  2  j 


(53) 


where  for  the  Expander  Types  1  and  4  of  course  Q  =  W  =  0 . 

Alternatively,  the  above  equations  can  be  replaced  with  the  following  equations 

•  Total  Mass  Balance: 


4U2  [(1-  e2 )  pG>2  +  e2  pf2  ]  =  mGX  +  m 


/.i 


(54) 


where: 


thG  X  ^V,  1  1  (^-0 

mG  =  mass  flow  rate  vapor  and  noncondensable 

mv  =  vapor  mass  flow  rate 

mn  =  noncondensable  mass  flow  rate 
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mf  - 


liquid  mass  flow  rate 


Note  that  everywhere: 


m„  =  mrx„ 

n  on 


Noncondensable  Mass  Balance: 


A7U2  [(1-  e2 )  pG2xn>2  ]  =  ihGxxnl  =  mnl  (57) 


Mixture  Energy  Balances: 


A^J2  (1  s2)/?G2  hG2  +  ^  U2  ^ + s2  Pj-  2  ^/,2  +  2  ^2 
=  Q~W  +  my  ihgl  +  mnlhn  l  +  mG]  hjx+  mflhfA 


Knowing  all  inlet  parameters,  and  assuming  that  P2  is  also  known,  Eqn.  (50-52),  or 
alternatively  Eqns.  (53-57),  with  appropriate  constitutive  relations,  must  be  solved  in 
order  to  obtain,  e2,hG  2,hf  2,U2,and  xn2 .  The  constitutive  relations  will  be  described  in 
the  forthcoming  Sect  4.2.3. 

For  the  Expander  Type  1,  at  inlet,  droplets  with  a  log-normal  size  distribution,  at  a 
total  mass  flow  rate  of  mf  l  and  at  a  total  number  rate  of  mf  i  and  at  a  total  number  rate 

of  Nl  X  are  injected  into  the  flow  field. 

The  log-normal  size  distribution,  which  is  likely  to  apply  to  spray  droplets,  can  be 
represented  as: 


N(D) 

N. 


1  (lnZ)-lnZ))2 

=  F(D)  = - ■F=^exp  - T-— - 

\n(jgyf2aD  2(ln<rg) 


(57-a) 


where: 
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N(D)dD  =  number  density  of  droplets  with  diameters  in  the  D  to  D  +  dD 


range 

Nt  =  total  number  density  of  droplets 

crg  =  geometric  standard  deviation 
obviously, 

oo 

e=  [^D’f(D)dD  (59) 

dL6 

The  droplet  (number)  average  and  volume  (or  mass)  average  diameter,  D ,  and  Dv ,  respectively, 
are  defined  as 

oo 

D  =  \f(D)DdD  (60) 

o 

Dy  =  J/  (D)D3dD  (61) 

o 

It  can  also  be  shown  that  for  the  log-normal  size  distribution: 
lnA,=ln.D  +  3(lncrg)2  (62) 

We  assume  that,  as  a  result  of  condensation,  the  droplets  grow  in  size  such  that  the  size 
distribution  function  does  not  change  (i.e.,  it  remains  log-normal,  with  the  geometric  standard 
deviation  unchanged).  As  a  result,  the  volume-average  droplet  diameter  at  exit  will  be 

(  fi  V/3 

BV1  =  -53-  (63) 

.  {*Nt,2  J 

Note  that  the  conservation  of  droplet  number  requires  that 
AxUxNt>x  =  A2U2Nt2  (64) 

An  equation  similar  to  Eqn.  (63)  can  obviously  be  also  written  for  the  inlet  to  the  expander. 


For  the  Expander  Type  4,  Eqns.  (53,-57)  apply,  except  that  now  mf  x  =  0 .  Also,  for  this 
expander  type,  the  droplet  volume-averaged  diameter  should  be  provided  by  the  user.  Knowing 
DV2,  the  number  density  and  number  flow  rate  can  be  easily  calculated. 

4.2.2  Expander  Types  2  and  3 

These  expander  types,  as  mentioned  earlier,  are  actually  adiabatic  turbine  expanders  with 
known  isentropic  efficiency  rjT ,  where: 

tjt  =  W/Ws  (65) 

where  W  and  Ws  represent  the  turbine  power,  and  its  isentropic  power,  respectively. 

For  these  expanders,  first  Ws  is  calculated.  The  calculated  W  is  then  used  in  Eqn.  (64), 

and  W  is  found.  Subsequently,  Egns.  (53-57),  along  with  the  constitutive  relations  describe  in 
the  next  section,  are  solved. 

To  obtain  Ws ,  the  following  is  done.  In  Eqn.  (57),  W  is  replaced  with  Ws .  Equations 
(53-57)  are  then  solved,  along  with  the  following  equation  which  ensures  an  isentropic  expansion 
of  the  fluid  mixture  in  the  expander: 

=  4^2  [(1-  e2K>2+  e2  s/>2]  (66) 

sfx ,  sg  l  =  saturated  water  and  vapor  entropies  at  inlet 

sf  2  =  saturated  water  entropy  at  exit  (for  the  vapor  partial  pressure  at  exit) 

sG  2  =  vapor-non  condensable  mixture  entropy 

Note  that 
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(67) 


SG,2  —  Xn,  2Sn,2  +  0  Xn,l)Sg,2 

Iterative  solution  of  the  above  equations,  along  with  other  appropriate  constitutive 
relations,  leads  to  the  calculation  of  Ws ,  as  well  as  the  isentropic  values  of  e2 ,  hG  2 ,  hf2 ,  etc. 
Only  Ws  is  useful,  however. 


4.2.3  Constitutive  Relations 

These  relations  apply  to  the  expander  inlet,  outlet,  or  anywhere  else  in  the  flow  field. 


P  =  Pv+Pn  (68) 

where 

P  =  Total  local  pressure 

PV,P„  =  Vapor  and  noncondensabe  partial  pressures,  respectively 


Pn  =  Ynp 

xJM_ „ 

"  Xn  ,  l~Xn 

Mn  My 


(70) 


(69) 


where: 


Y  =  mole  fraction  of  noncondensable 

ft 

Mn ,  Mv  =  molar  masses  of  noncondensable  and  vapor  respectively 


T  =  Tsat{Py)  (71) 

A,-P,W+-3*-  (72) 

—  T 
Mn 

R  =  universal  gas  constant 
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p  =  saturated  vapor  density 
hG=(l-xn)hg+xnh„(T)  (73) 


where: 


hn  =  Specific  enthalpy  of  the  noncondensable  species 
Also,  since  we  treat  the  noncondensabe  (air  in  this  care)  as  an  ideal  gas,  can  write: 


(74) 


where s°(T2)  and  s°n{Tx)  are  the  following  standard  function 


<(n=jc„^  (75) 


Values  of  s°(T)  for  air  were  calculated  by  interpolation  using  the  tables  of  Black  and  Hartley 
(1996). 


4.2.4  Solution  Method 

For  each  expander  type,  the  relevant  set  of  algebraic  equation  for  is  a  closed,  non-linear 
equation  system.  There  equations  are  solved  using  the  Newton  iteration  method.  Briefly,  for  the 
equation: 

f(x)  =  0  (76) 

the  iterations  in  search  of  the  roots  of  Eqn.  (75)follow: 

*«+i  =*«-/(*.  V/W  (77) 
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where  subscripts  n  and  n  + 1  represent  their  respective  iteration  levels,  and  /'  is  the  derivative 
of  / (x)  with  respect  to  x . 

When  a  set  of  nonlinear  equations  are  dealt  with,  Eqn.(76)  is  replaced  with 
Xn+1=Xn-F(Xn)/J  (78) 

where  X  now  represents  a  column  vector  with  elements  xvx2...,and  F{X)  is  a  column  vector 

with  elements  /(x,,x2, ....), /2(x„x2, ....),  etc.,  and  parameter  J  is  the  determined  of  the  Jacobian 

matrix.  The  elements  of  the  Jacobian  matrix  are  themselves  approximately  obtained  from 

df.  fiixj+m-ftixj) 
dxj  H 

where  H  is  a  small  number. 
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Fig.  8  Schematic  of  a  cyclone  with  Tangential  Gas  Inlet  (Leith  and  Licht,  1972) 
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4.3  Cyclone-Type  Separator 
4.3.1  General  Remarks 

Cyclone  separators  are  devices  that  partially  remove  particles  suspended  in  a  particle¬ 
laden  gas  or  liquid.  Particle  removal  is  achieved  by  providing  a  highly  curved  flow  path  for  the 
fluid  mixture.  Unless  the  particles  are  neutrally-buoyant,  they  would  deviate  from  the  carrier 
fluid’s  stream  lines  due  to  their  inertia,  and  particle-wall  impactions  would  result  in  the  removal 
of  particles  from  the  mixture. 

Cyclone  separators  and  their  hydrodynamic  and  particle  transport  characteristics  have 
been  extensively  studied  in  the  past.  Earlier  modelers  attempted  to  derive  closed  mathematical 
expressions  or  simple  numerical  methods  for  calculating  the  removal  or  collection  efficiency  of 
cyclones  with  known  geometric  configurations  CFD-type  models,  which  can  address  cyclone 
separators  with  rather  arbitrary  geometric  configurations,  have  been  developed  and  published 
more  recently.  Virtually  all  published  studies,  however,  are  concerned  with  particles  that  do  not 
undergo  phase  change.  A  brief  review  of  the  recent  literature  follows. 

An  important,  widely-used  simple  model  for  the  collection  efficiency  of  cyclones  is  the 
model  of  Leith  and  Licht  (1972).  The  configuration  of  a  cyclone  with  tangential  gas  inlet, 
modeled  by  the  latter  authors,  is  shown  schematically  in  Fig.  8.  The  semi-mechanistic  model  of 
Leith  and  Licht  can  also  be  applied  to  cyclones  with  circular  inlets.  Other,  more  recent  studies 
dealing  with  cyclones  include  the  following.  Comas  et  al.  (1991)  experimentally  studied  the 
effect  of  dust  and  inlet  vanes  on  cyclone  collection  efficiency  and  pressure  drop.  Dust  particles 
reduced  pressure  drop,  while  vanes  reduced  the  removal  efficiency  and  pressure  drop  both. 
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Dyakowski  and  Williams  (1993)  numerically  modeled  the  flow  in  narrow  hydrocyclones, 
using  the  k-  e  turbulence  model.  Griffith  and  Boysan  (1996)  argued  that  cyclone  models  such 
as  the  model  of  Leith  and  Licht  (1972)  are  good  for  large,  industrial  cyclones,  and  may  not  be 
appropriate  for  small  units.  They  performed  CFD-type  analysis  of  three  model  cyclone  types 
and  assessed  the  performance  of  three  relevant  empirical  models.  Hoekstra  et  al.  (1999) 
performed  an  experimental  and  numerical  study  of  turbulent  flow  in  a  gas  cyclone.  In  their  CFD 
analysis,  they  compared  the  performance  of  three  different  turbulence  models,  showing  that  the 
RSM  performed  better  than  the  k-  e  and  RNG  models.  Chmielniak  and  Bryczkowski  (2000) 
presented  a  theoretical  model  for  a  patented  cyclone  device  which  employs  swirling  baffles  and  a 
bottom  clean  gas  take-off.  Their  theoretical  model  was  based  on  the  modification  of  the 
aforementioned  model  of  Leicht  and  Leith  (1972). 

A  simple,  one-dimensional  model  based  on  the  assumption  of  a  spiral  motion  with 
constant  acceleration  in  the  cyclone  has  also  been  suggested  recently  (Avci  and  Karagoz,  2000). 
However,  it  has  not  been  sufficiently  validated  against  experimental  data. 

As  noted  earlier,  the  above  studies  as  well  as  others  in  the  open  literature,  are  primarily 
concerned  with  particles  without  phase  change.  The  physical  situation  we  deal  with  is  different, 
however.  The  vapor-noncondensable-droplet  mixture  undergoes  depressurization  in  the  cyclone, 
which  results  in  further  condensation  of  vapor  on  droplets,  the  growth  of  droplets  size,  and  the 
reduction  of  the  mixture  mean  velocity.  In  addition,  the  geometric  configuration  of  the  cyclone 
of  interest  may  not  exactly  coincide  with  the  standard  cyclone  depicted  in  Fig.  8. 

Detailed  CFD-type  analysis  is  evidently  needed  for  cases  of  interest  here,  and  should  be 
performed  in  the  future.  At  this  point,  however,  in  light  of  the  uncertainties  regarding  the  final 
separator  design,  only  a  simple  model  will  be  developed. 
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4.3.2  Cyclone  Model 

The  cyclone  separator  model  will  be  based  on  modification  of  the  model  of  Leith  and 
Licht  (1972),  in  order  to  account  for:  (a)  the  occurrence  of  condensation  during  the  flow  in  the 
cyclone;  and  (b)  the  geometric  difference  (primarily  the  absence  of  the  comic  part,  see  Fig.  8.) 
between  the  currently  selected  design  and  the  cyclone  model  by  the  latter  authors. 

The  modified  model,  in  addition  to  the  assumptions  that  are  the  bases  of  the  model  of 
Leith  and  Licht  (1972),  assumes  the  following. 

1.  At  inlet  to  the  cyclone  the  particle  size  distribution,  along  with  other  flow 
parameters,  are  all  known. 

2.  The  vapor-noncondensable-droplet  mixture  is  at  thermodynamic  equilibrium 
everywhere  in  the  system. 

3.  For  particle  removal  calculations  the  particles  can  be  treated  as  monodisperse, 
with  the  volume-average  diameter  as  their  effective  diameter. 

4.  The  removal  (collection)  efficiency  of  the  cyclone  can  be  obtained  using  the 
method  of  Leith  and  Licht,  based  on  the  particle  size  characteristics  midway  in  the 
cyclone. 

5.  The  collected  droplets  all  leave  the  cyclone  as  saturated  liquid  corresponding  to 
the  cyclone  exit  pressure. 

The  cyclone  performance  is  accordingly  obtained  by  the  simultaneous  solution  of 
overall  mass  and  energy  balance  equation,  along  with  the  expressions  from  Leith  and  Licht 
model.  The  set  of  coupled  equations  are  as  follows. 

•  Mixture  Mass  Balance 
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(80) 


Km  +  AU2  [(1-  e2  )  Pc,  2+  e2  Pf,l\=  <,]  +  '”/,! 

•  Noncondensable  Mass  Balance 

A2U2[(l-s2)pG2x„2]  =  mnA 

•  Mixture  Energy  Balance 


™remhf,  2+4^2 


Pc,2  (l-e2) 


+  e2  pf,2(hfa  +|[/22 )]  =  <A.i  +  KK\ 

+Ai  "2  ^ 


(81) 


(82) 


where  the  removed  liquid  flow  rate  is  given  by: 


mrm=(N,AUl-NtU,A1)pf-D, 


ft  ^3 
V 


(83) 


where  iV,  and  tV2  are  droplet  number  densities  at  inlet  and  exit,  respectively.  Also, 


A  ~  2  {Pvi +  DVs ) 


(84) 


(85) 


6,=  *,^. 


(86) 


The  collection  efficiency,  77,  is  now  formulated  following  Leith  and  Licht  (1972) 
(see  Fig.8  for  the  definition  of  geometric  parameters  used): 


ij  =  1  -  exp 


-2  c(y/) 


1 

2n+2 


KK 


(87) 


(88) 
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(89) 


¥  = 


18  piD 


(«+l) 


where  Vs  is  the  annular  shape  volume  above  the  exit  duct  through  the  mid-level  of  the 
entrance  duct,  and: 


11 

(90) 

V„=(h-s)~D> 

(91) 

(12  XL)0'14 
n  = - - - 

(92) 

2.5 
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k=± 
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H 
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It  should  be  mentioned  that  in  the  method  of  Leith  and  Licht,  a  natural  length,  1 ,  and  a 
corresponding  volume  Vnl  and  kc  are  defined  as: 


/  =  23De(D2  /  ab)m 

(97) 

V  i  =  l—D2 

4 

(98) 

K  = 


V.+  —  V. 
s_  2 

D3 


(99) 
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Since  our  current  separator  design  does  not  included  a  conical  segment,  we  have  defined 
kc  according  to  Eqn.  (94). 

Also,  when  the  inlet  to  the  cyclone  is  circular  rather  than  rectangular,  the  above 
equations  apply  as  long  as  ab  is  replaced  with  the  cross-sectional  area  of  the  inlet 

channel,  and  kakb  is  replaced  with  n  / 4  (D;  ID)1 ,  with  D}  representing  the  inlet  channel 


diameter. 
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Section  1 
Introduction 


1.1  Problem  Background 

EVAPCM  is  a  computer  code  that  has  been  developed  to  mechanistically  model  the 
evaporator  process  in  a  uniform  cross-section  channel.  EVAPMC  is  based  on  the 
numerical  solution  of  sets  of  one-dimensional  conservation  equations  in  various  segments 
of  a  boiling  channel  (the  subcooled  segment  near  the  inlet,  the  two-phase  boiling 
segment,  etc.),  everywhere  using  appropriate  closure  relations,  and  is  capable  of 
providing  hydrodynamic  details  of  the  flow  field  it  models. 

EVAPMCcontains  built-in  closure  relations  that  are  suitable  for  the  analysis  of  the 
evaporation  of  water  in  a  small  flow  passage  with  rectangular  cross-section.  Extensive 
property  routines  for  water  are  included  in  the  code.  The  mechanistic  modeling 
approach,  however,  is  quite  general,  and  constitutive  relations  suitable  for  other 
configurations  and  working  conditions  can  be  easily  incorporated  into  the  code. 


1.2  Report  Summary 

Section  2  briefly  describes  the  phenomena  that  are  modeled  in  this  computer  code.  The 
code  structure  along  with  its  flow  diagram  is  presented  in  Section  3.  Section  4  briefly 
describes  each  of  the  routines  used  in  the  code.  The  input  and  output  are  described  in 
Section  5.  A  sample  case  is  provided  in  Section  6,  where  a  complete  listing  of  all  the 
input  and  output  is  included.  The  references  are  listed  in  Section  7. 
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Section  2 

Description  of  the  Theory 

The  modeled  flow  channels  can  be  divided  into  4  segments,  as  shown  in  the  schematic  in 
Figure  2.1.  The  model  uses  a  separate  set  of  steady-state,  one-dimensional  conservation 
equations  in  each  region. 

In  region  I,  the  subcooled  liquid  region,  momentum  and  energy  conservation  ordinary 
differential  equations  are  numerically  solved,  using  pressure  and  enthalpy  as  the  state 
variables,  until  saturation  is  reached,  i.e.,  until  the  fluid  enthalpy  falls  below  the  saturated 
liquid  enthalpy  associated  with  the  local  pressure.  To  avoid  numerical  problems, 
however,  we  allow  the  code  to  slightly  “undershoot”.  Thus,  evaporation  is  assumed  to 
start  when  the  local  equilibrium  vapor  quality  corresponds  to  a  homogeneous  void 
fraction  of  0.25%.  Once  the  latter  condition  is  satisfied,  Region  II  starts. 

Region  II  represents  the  segment  of  the  channel  where  two-phase  flow  with  various 
familiar  regimes  takes  place.  In  this  region  the  separated  (two-fluid)  conservation 
equations  are  numerically  solved.  It  is  assumed  that  the  liquid  and  vapor  phases  remain 
at  thermodynamic  equilibrium,  while  maintaining  different  velocities.  A  simple  flow 
regime  map,  composed  of  bubbly,  intermittent  (chum-slug)  and  annular  flow  patterns,  is 
used  for  the  calculation  of  the  interfacial  forces.  The  parameters  in  each  flow  regime 
(e.g.,  bubble  size  in  bubbly  flow)  are  currently  represented  with  the  rectangular  channel 
restrictions  in  mind.  These  parameters,  however,  can  be  easily  modified.  The  ordinary 
differential  conservation  equations  are  numerically  integrated,  with  pressure,  vapor 
velocity,  liquid  velocity,  and  void  fraction,  as  the  state  variables.  Region  II  is  assumed  to 
end  when  the  local  equilibrium  quality  exceeds  20%,  and  the  local  void  fraction  predicted 
by  the  two-fluid  model  exceeds  95%. 

Region  III  represents  evaporation  in  a  high-quality  two-phase  flow  pattern,  which,  based 
on  the  phenomenology  of  boiling  in  large  channels,  is  likely  to  represent  the  annular  or 
annular/mist  flow  pattern.  In  this  region  thermodynamic  equilibrium  is  assumed  between 
the  two  phases,  and  the  mixture  mass,  momentum,  and  energy  ordinary  differential 
equations  are  numerically  solved,  with  pressure,  equilibrium  quality,  and  vapor  velocity. 
The  velocity  slip  between  the  two  phases  is  accounted  for  by  using  a  correlation  due  to 
Zivi  (1964),  which  has  been  derived  based  on  the  assumption  of  minimum  entropy 
generation  in  annular  flow.  Region  III  is  assumed  to  end  when  the  local  equilibrium 
quality  reaches  0.999. 

Region  IV  represents  superheated  pure  vapor.  The  conservation  equations  and  the 
solution  method  are  similar  to  Region  I,  with  compressibility  of  the  vapor  phase  being  the 
major  difference  between  the  two  cases. 

It  should  be  mentioned  that  not  all  the  aforementioned  regions  have  to  occur  in  a 
simulation  run.  Thus,  for  example,  Region  IV  will  not  occur  if  the  exit  of  the  channel  is 
reached  before  the  equilibrium  quality  reaches  0.999. 
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Figure  2.1  Schematic  of  the  Modeled  Channel. 


Section  3 
Code  Structure 


The  logical  flow  of  the  computer  code  is  shown  schematically  in  Figure  3-1 .  The 
program  begins  by  reading  and  writing  the  input  data.  Following  this,  the  variables  that 
depend  on  the  axial  coordinate,  z,  are  initialized.  Then  the  program  begins  marching 
down  the  axial  coordinate,  in  the  direction  of  increasing  z.  Marching  along  z  in  the 
aforementioned  Region  I  (See  Figure  2.1)  is  performed  in  subroutine  MARCH.  Once 
the  criterion  for  the  onset  of  two-phase  flow  is  reached,  the  program  switches  to 
subroutine  MARCH2,  where  Regions  II  and  III  are  addressed.  Finally,  once  complete 
evaporation  occurs,  the  program  switches  to  subroutine  MARCH3,  where  Region  IV  is 
addressed,  until  the  channel  exit  is  reached.. 

A  summary  of  the  calculation  results  is  printed  at  user-specified  intervals. 
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Section  4 

Description  of  Subroutines  and  Functions 

This  section  contains  a  description  of  each  subroutine  and  function  contained  in  this 
computer  code. 


Subroutines 


CPLT 

DVODE 

FRIEDEL 

FUNC 

FUNC2 

FUNC2F 

FUNC3 

INVER 

MARCH 

MARCH2 

MARCH3 

POLINT 

PRPSAT 


Finds  the  subcooled  liquid  water  pecific  heat  given  the  pressure  and 
temperature. 

Solves  initial-value  problems  represented  by  sets  of  ordinary  differential 
equations.  Stiff  and  non-stiff  systems  can  both  be  efficiently  solved. 

Calculates  the  two-phase  wall  friction  using  Friedel’s  (1977)  model.  (This 
subroutine  is  not  used  at  this  time.) 

Calculates  and  sets  up  the  matrix  elements,  and  performs  necessary 
matrix  inversions,  for  the  ode’s  to  be  solved  in  MARCH  subroutine. 

Calculates  and  sets  up  the  matrix  elements,  and  performs  necessary  matrix 
inversions,  for  the  ode’s  to  be  solved  in  MARCH2  subroutine  for  Region 
II. 

Calculates  and  sets  up  the  matrix  elements,  and  performs  necessary  matrix 
inversions,  for  the  ode’s  to  be  solved  in  MARCH2  subroutine  for  Region 

m. 

Calculates  and  sets  up  the  matrix  elements,  and  performs  necessary  matrix 
inversions,  for  the  ode’s  to  be  solved  in  MARCH3  subroutine. 

This  subroutine  inverts  a  matrix  using  a  standard  Gauss-Jordan 
elimination  technique. 

The  main  routine  for  Region  I  calculations. 

Main  routine  for  the  calculations  for  Regions  II  and  III. 

Main  routine  for  the  calculations  for  Region  IV. 

This  subroutine  is  a  Lagrangian  interpolation  routine,  and  is  used  to 
interpolate  between  thermodynamic  values  in  property  tables. 

A  block  of  subroutines  and  functions  that  calculate  various  water  and 
steam  properties  at  saturated  conditions. 
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PSATF 

RKAMS 


THCND 

THCONV 

VISCOL 

VISCOV 


Finds  the  saturation  pressure  given  the  temperature. 

This  subroutine  is  a  fixed-step,  4th  order  Runge-Kutta  solver  for  ordinary 
differential  equations.  It  is  utilized  in  MARCH,  MARCH2  and  MARCH3 
routines  for  the  integration  of  all  of  the  ordinary  differential  equation  sets. 
(Not  used.) 

Finds  the  thermal  conductivity  of  water  using  temperature  and  density. 
Finds  vapor  thermal  conductivity  given  pressure  and  temperature. 

Finds  the  viscosity  of  water  given  the  saturation  temperature. 

Finds  vapor  viscosity  given  pressure  and  temperature. 


Functions 

CPA28 

CPGAS 

CPVPT 

DHGPSAT 

EKNTHE 

EMUNTHE 

HAT28 

HENRY 

HFPSAT 

HGPSAT 

SIGTF 

TSATF 

YFT 

VGPT 


Finds  the  specific  heat  of  N2  gas  (not  used  here). 

Finds  the  specific  heat  at  constant  pressure  for  superheated  steam  given 
the  temperature. 

Finds  the  specific  heat  of  water  vapor  given  pressure  and  temperature. 
Finds  the  derivative  of  enthalpy  versus  pressure  for  water  vapor. 

Finds  thermal  conductivity  for  non-condensible  He,  given  the  temperature. 
Finds  viscosity  for  non-condensible  He,  given  the  temperature. 

Finds  enthalpy  for  non-condensible  N2  given  the  temperature. 

Finds  Henry’s  constant  for  N2  and  water.  (This  function  is  not  used.) 

Finds  the  liquid  saturation  enthalpy  given  the  pressure. 

Finds  the  vapor  saturation  enthalpy  given  the  pressure. 

Finds  the  surface  tension  of  water. 

Given  the  pressure,  this  returns  the  saturation  temperature. 

Given  the  saturation  temperature,  this  function  finds  the  density  of  water. 
Finds  specific  volume  given  pressure  and  temperature. 


8 


Section  5 

Input  and  Output 

This  section  describes  the  input  parameters  necessary  to  run  this  computer  code  and  an 
explanation  of  all  the  associated  output. 

INPUT 

Table  5.1  defines  all  the  input  parameters.  The  inputs  are  divided  into  five  basic  groups, 
as  follows: 

1 .  Geometric  parameters 

2.  Inlet  conditions 

3.  Heat  input 

4.  Pressure  drop  parameters 

5.  Integration  and  Printing  parameters 

Each  input  variable  is  on  a  separate  line  in  the  input  file.  The  variables  must  be  input  in 
the  following  order. 

Line  1:  MN,  MV,  WIDE,  HIGH,  TUBEL,  THETA 

Line  2:  TIN,PIN,VIN 

Line  3:  QW,  ZHEATS,  ZHEATF 

Line  4:  ELOSSIRR,  ELOSSREV 

Line  5:  BETA 

Line  6:  ZPRNT 
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Table  5.1 
Input  Definition 

1.  Geometric  Parameters 


Parameter 

Units 

Description 

HIGH 

Meters 

Height  of  channel  cross  section 

WIDE 

Meters 

Width  of  channel  cross-section 

THETA 

Degrees 

Angle  of  orientation  from  horizontal 

TUBEL 

Meters 

Channel  length 

2.  Inlet  Conditions 

Parameter 

Units 

Description 

BETA 

1/Degre  K 

Liquid  thermal  expansion  coefficient 

MN 

kg/kgmole 

Molecular  weight  of  noncondensibles 

MV 

kg/kgmole 

Molecular  weight  of  water  vapor 

PIN 

Pa 

Inlet  pressure  of  fluid 

TIN 

°C 

Inlet  temperature  of  fluid 

YIN 

m/s 

Inlet  velocity  of  fluid 

3.  Heat  Input 

Parameter 

Units 

Description 

QW 

W/m 

Heat  input  rate  per  unit  of  channel  length 

ZHEATS 

M 

Location  along  channel  where  heat  input 
begins 

ZHEATF 

M 

Location  along  channel  where  heat  input 

ends 
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4.  Pressure  Drop  Parameters 


Parameter _ Units _ Description _ 

ELOSSIRR  -  Irreversible  Loss  Coefficient  at  Inlet 


ELOSSREV  -  Reversible  Pressure  drop  Coefficient  at 

Inlet 


5.  Integration  and  Printing  Parameters 

Parameter _ Units _ Description 

ZPRNT  m  Printing  Increment 
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Output  Information 
Includes: 


1 .  All  input  variables 

2.  Local  conditions  along  length  of  channel,  divided  by  flow  regimes.  Pressure  and 
enthalpy  are  given  at  all  locations.  At  the  locations  where  two  phase  flow  occurs, 
the  quality,  x,  and  void  fraction,  a,  are  given  as  well  as  the  gas  and  liquid 
velocities,  ug  and  Uf,  respectively,  where  appropriate. 

3.  Location  of  onset  of  boiling,  switch  to  film  flow  calculations  and  location  of 
single-phase  vapor  flow  are  printed  out. 
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Section  6 

Sample  Input  and  Output 


Sample  Input 

29.0,18.0,0.005,0.06,0.21,0.0, 

20.0,2. 03e6, 0.3, 

1.5E6,.01,0.2, 

.3, .5, 

304.E-6, 

0.001 

Sample  Output 

PROGRAM  INPUT  DATA 
Molecular  Weights  of  Noncondensibles  and  Water 

MN  =  29.0000  MV  =  18.0000 

Channel  Width  Channel  Height  Length  Angle  of  Inclination 
.0050  .0600  .2100  .0000 

Inlet  Conditions 

Temperature  Pressure  Velocity 
(Deg  C)  (Pa)  (m/s) 

20.000000  2030000.000000  .300000 


Heat  Input  =  .150000E+07  W/m 
Heating  starts  at  z  =  .010000  m 
Heating  ends  at  z  =  .200000  m 

Inlet  irreversible  loss  coefficient  =  .300 
Inlet  reversible  loss  coefficient  =  .500 
Thermal  expansion  coef.  of  liquid  =  3.04000E-04  1/K 
Print  interval  =  1.00000E-03  m 


Liquid  Single  Phase  Results 

z  (m)  T  (K)  P  (Pa)  Xeq  (-) 


.00000  293.15 
.00100  321.95 
.00200  321.95 
.00300  321.95 


2.0300E+06  -3.0963E-01 
2.0300E+06  -3.9534E-01 
2.0300E+06  -3.9534E-01 
2.0300E+06  -3.9534E-01 
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.00400  321.95 
.00500  321.95 
.00600  321.95 
.00700  321.95 
.00800  321.95 
.00900  321.95 
.01000  321.95 
.01100  325.62 
.01200  329.30 
.01300  332.98 
.01400  336.65 
.01500  340.33 
.01600  344.00 
.01700  347.68 
.01800  351.36 
.01900  355.03 
.02000  358.71 
.02100  362.38 
.02200  366.06 
.02300  369.74 
.02400  373.41 
.02500  377.09 
.02600  380.76 
.02700  384.44 
.02800  388.12 
.02900  391.79 
.03000  395.47 
.03100  399.14 
.03200  402.82 
.03300  406.50 
.03400  410.17 
.03500  413.85 
.03600  417.52 
.03700  421.20 
.03800  424.88 
.03900  428.55 
.04000  432.23 
.04100  435.90 
.04200  439.58 
.04300  443.26 
.04400  446.93 
.04500  450.61 
.04600  454.28 
.04700  457.96 
.04800  461.64 
.04900  465.31 


2.0300E+06  -3.9534E-01 
2.0300E+06  -3.9534E-01 
2.0300E+06  -3.9534E-01 
2.0300E+06  -3.9534E-01 
2.0300E+06  -3.9534E-01 
2.0300E+06  -3.9534E-01 
2.0300E+06  -3.9534E-01 
2.0300E+06  -3.8651E-01 
2.0300E+06  -3.7768E-01 
2.0300E+06  -3.6885E-01 
2.0300E+06  -3.6002E-01 
2.0300E+06  -3.51 19E-01 
2.0300E+06  -3.4237E-01 
2.0300E+06  -3.3354E-01 
2.0300E+06  -3.2471E-01 
2.0300E+06  -3.1588E-01 
2.0300E+06  -3.0705E-01 
2.0300E+06  -2.9822E-01 
2.0300E+06  -2.8939E-01 
2.0300E+06  -2.8056E-01 
2.0300E+06  -2.7173E-01 
2.0300E+06  -2.6290E-01 
2.0300E+06  -2.5407E-01 
2.0300E+06  -2.4524E-01 
2.0300E+06  -2.364 IE-01 
2.0300E+06  -2.2758E-01 
2.0300E+06  -2.1875E-01 
2.0300E+06  -2.0993E-01 
2.0300E+06  -2.0110E-01 
2.0300E+06  -1.9227E-01 
2.0300E+06  -1.8344E-01 
2.0300E+06  -1.7461E-01 
2.0300E+06  -1.6578E-01 
2.0300E+06  -1.5695E-01 
2.0300E+06  -1.4812E-01 
2.0300E+06  -1.3929E-01 
2.0300E+06  -1.3046E-01 
2.0300E+06  -1.2163E-01 
2.0300E+06  -1.1280E-01 
2.0300E+06  -1.0397E-01 
2.0300E+06  -9.5144E-02 
2.0299E+06  -8.6314E-02 
2.0299E+06  -7.7485E-02 
2.0299E+06  -6.8655E-02 
2.0299E+06  -5.9826E-02 
2.0299E+06  -5.0997E-02 
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.05000  468.99 
.05100  472.66 
.05200  476.34 
.05300  480.02 
.05400  483.69 


2.0299E+06  -4.2167E-02 
2.0299E+06  -3.3338E-02 
2.0299E+06  -2.4509E-02 
2.0299E+06  -1.5679E-02 
2.0299E+06  -6.8500E-03 


**  Boiling  Started  at  z  =  .0550  (m) 
Boiling  Region  Results 


z(m) 

P  (Pa)  Xeq  (-)  Void  (-)  Ug  (m/s) 

Uf  (m/s) 

.05600 

2.0299E+06 

1.0923E-02 

4.1272E-01 

7.8134E-01 

.05700 

2.0298E+06 

1.9814E-02 

5.3857E-01 

1.0880E+00 

.05800 

2.0298E+06 

2.8621E-02 

6.1971E-01 

1.3707E+00 

.05900 

2.0297E+06 

3.7630E-02 

6.7914E-01 

1.6387E+00 

.06000 

2.0297E+06 

4.6165E-02 

7.2427E-01 

1.9004E+00 

.06100 

2.0296E+06 

5.5174E-02 

7.5843E-01 

2.1623E+00 

.06200 

2.0296E+06 

6.4183E-02 

7.8357E-01 

2.4292E+00 

.06300 

2.0296E+06 

7.2704E-02 

8.0019E-01 

2.7081E+00 

.06400 

2.0295E+06 

8.1729E-02 

8.1164E-01 

2.9946E+00 

.06500 

2.0295E+06 

9.1415E-02 

8.2191E-01 

3.2778E+00 

.06600 

2.0295E+06 

9.9835E-02 

8.3116E-01 

3.5584E+00 

.06700 

2.0294E+06 

1.0909E-01 

8.3952E-01 

3.8369E+00 

.06800 

2.0294E+06 

1.1834E-01 

8.471  IE-01 

4.1137E+00 

.06900 

2.0294E+06 

1.2574E-01 

8.5403E-01 

4.3890E+00 

.07000 

2.0293E+06 

1.3499E-01 

8.6036E-01 

4.6630E+00 

.07100 

2.0293E+06 

1.4424E-01 

8.6617E-01 

4.9360E+00 

.07200 

2.0292E+06 

1.5259E-01 

8.7153E-01 

5.2082E+00 

.07300 

2.0292E+06 

1.6375E-01 

8.7648E-01 

5.4795E+00 

.07400 

2.0292E+06 

1.7212E-01 

8.8108E-01 

5.7501E+00 

.07500 

2.0291E+06 

1.8050E-01 

8.8535E-01 

6.0201  E+00 

.07600 

2.029 1E+06 

1.8887E-01 

8.8933E-01 

6.2896E+00 

.07700 

2.0290E+06 

1.9724E-01 

8.9306E-01 

6.5586E+00 

.07800 

2.0290E+06 

2.0561E-01 

8.9655E-01 

6. 8271 E+00 

.07900 

2.0290E+06 

2.1678E-01 

8.9982E-01 

7.0953E+00 

.08000 

2.0289E+06 

2.2515E-01 

9.029 IE-01 

7.3631  E+00 

.08100 

2.0289E+06 

2.3359E-01 

9.0581E-01 

7.6306E+00 

.08200 

2.0288E+06 

2.4203E-01 

9.0856E-01 

7.8978E+00 

.08300 

2.0288E+06 

2.5047E-01 

9.1115E-01 

8.1647E+00 

.08400 

2.0287E+06 

2.5891E-01 

9.1361E-01 

8.4314E+00 

.08500 

2.0287E+06 

2.6735E-01 

9.1594E-01 

8.6979E+00 

.08600 

2.0287E+06 

2.8001E-01 

9.1815E-01 

8.9642E+00 

.08700 

2.0286E+06 

2.8845E-01 

9.2026E-01 

9.2303E+00 

.08800 

2.0286E+06 

2.9689E-01 

9.2226E-01 

9.4963E+00 

.08900 

2.0285E+06 

3.0532E-01 

9.2417E-01 

9.7621E+00 

.09000 

2.0285E+06 

3.1376E-01 

9.2599E-01 

1.0028E+01 

.09100 

2.0284E+06 

3.2220E-01 

9.2773E-01 

1.0293E+01 

5.9404E-01 

7.493  IE-01 

9.0098E-01 

1.0582E+00 

1.2200E+00 

1.3796E+00 

1.5255E+00 

1.6368E+00 

1.7198E+00 

1.8015E+00 

1.8817E+00 

1.9603E+00 

2.0372E+00 

2.1124E+00 

2.1859E+00 

2.2576E+00 

2.3275E+00 

2.3956E+00 

2.4619E+00 

2.5265E+00 

2.5893E+00 

2.6503E+00 

2.7096E+00 

2.7671E+00 

2.8229E+00 

2.8769E+00 

2.9291E+00 

2.9796E+00 

3.0283E+00 

3.0752E+00 

3.1203E+00 

3.1636E+00 

3.205 1E+00 

3.2447E+00 

3.2824E+00 

3.3183E+00 
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.09200  2.0284E+06  3.3064E-01  9.2939E-01  1.0559E+01 

.09300  2.0283E+06  3.3908E-01  9.3098E-01  1.0824E+01 

.09400  2.0282E+06  3.4971E-01  9.3250E-01  1.1089E+01 

.09500  2.0282E+06  3.5612E-01  9.3396E-01  1.1355E+01 

.09600  2.028 1E+06  3.6895E-01  9.3535E-01  1.1620E+01 

.09700  2.028 1E+06  3.7536E-01  9.3669E-01  1.1885E+01 

.09800  2.0280E+06  3.8819E-01  9.3797E-01  1.2150E+01 

.09900  2.0280E+06  3.9460E-01  9.3920E-01  1.2415E+01 

.10000  2.0279E+06  4.0101E-01  9.4038E-01  1.2680E+01 

.10100  2.0279E+06  4.1384E-01  9.4151E-01  1.2946E+01 

.10200  2.0278E+06  4.2025E-01  9.4260E-01  1.321 1E+01 

.10300  2.0277E+06  4.2666E-01  9.4364E-01  1.3476E+01 

.10400  2.0277E+06  4.3948E-01  9.4464E-01  1.3741E+01 

.10500  2.0276E+06  4.4589E-01  9.4560E-01  1.4006E+01 

.10600  2.0275E+06  4.5872E-01  9.4652E-01  1.4272E+01 

.10700  2.0275E+06  4.6513E-01  9.4740E-01  1.4537E+01 

.10800  2.0274E+06  4.7154E-01  9.4825E-01  1.4803E+01 

.10900  2.0273E+06  4.8437E-01  9.4906E-01  1.5068E+01 

.11000  2.0273E+06  4.9078E-01  9.4983E-01  1.5334E+01 


3.3523E+00 
3.3843E+00 
3.4144E+00 
3.4425E+00 
3.4686E+00 
3.4927E+00 
3.5146E+00 
3.5345E+00 
3.5522E+00 
3.5677E+00 
3.5809E+00 
3.5919E+00 
3.6005E+00 
3.606 7E+00 
3.6105E+00 
3.6118E+00 
3.6106E+00 
3.6067E+00 
3.6002E+00 


Switch  to  F ilm  Flow  at  z  =  .111 0(m) 
Film  Region  Results 


z(m) 

P  (Pa)  Xeq  (-)  Void  (- 

-)  Ug  (m/s) 

Uf  (m/s) 

.11100 

2.027 1E+06 

5.0601E-01 

9.5223E-01 

1.5894E+01 

.11200 

2.0271E+06 

5.1484E-01 

9.5386E-01 

1.6171E+01 

.11300 

2.0270E+06 

5.2366E-01 

9.5527E-01 

1.6449E+01 

.11400 

2.0269E+06 

5.3249E-01 

9.573  IE-01 

1.6727E+01 

.11500 

2.0269E+06 

5.4132E-01 

9.5904E-01 

1.7005E+01 

.11600 

2.0268E+06 

5.5014E-01 

9.6008E-01 

1.7282E+01 

.11700 

2.0268E+06 

5.5897E-01 

9.6205E-01 

1.7560E+01 

.11800 

2.0267E+06 

5.6780E-01 

9.6298E-01 

1.7838E+01 

.11900 

2.0266E+06 

5.7662E-01 

9.6390E-01 

1.8116E+01 

.12000 

2.0265E+06 

5.8545E-01 

9.6514E-01 

1.8394E+01 

.12100 

2.0265E+06 

5.9427E-01 

9.6694E-01 

1.8671E+01 

.12200 

2.0264E+06 

6.0310E-01 

9.6788E-01 

1.8949E+01 

.12300 

2.0263E+06 

6.1192E-01 

9.6881E-01 

1.9227E+01 

.12400 

2.0263E+06 

6.2075E-01 

9.7060E-01 

1.9505E+01 

.12500 

2.0262E+06 

6.2958E-01 

9.7147E-01 

1.9783E+01 

.12600 

2.0261E+06 

6.3840E-01 

9.7231E-01 

2.0061E+01 

.12700 

2.0260E+06 

6.4723E-01 

9.7395E-01 

2.0339E+01 

.12800 

2.0260E+06 

6.5605E-01 

9.7475E-01 

2.0617E+01 

.12900 

2.0259E+06 

6.6488E-01 

9.7552E-01 

2.0895E+01 

.13000 

2.0258E+06 

6.7370E-01 

9.7675E-01 

2.1173E+01 

.13100 

2.0257E+06 

6.8253E-01 

9.7756E-01 

2.1451E+01 

3.6368E+00 
3.699 1E+00 
3.7547E+00 
3.8380E+00 
3.9114E+00 
3.9571E+00 
4.0466E+00 
4.0904E+00 
4.1343E+00 
4.1953E+00 
4.2870E+00 
4.3369E+00 
4.3869E+00 
4.4869E+00 
4.5369E+00 
4.5869E+00 
4.6868E+00 
4.7368E+00 
4.7868E+00 
4.8676E+00 
4.9228E+00 
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.13200  2.0256E+06  6.9135E-01 
.13300  2.0256E+06  7.0017E-01 
.13400  2.0255E+06  7.0900E-01 
.13500  2.0254E+06  7.1782E-01 
.13600  2.0253E+06  7.2665E-01 
.13700  2.0252E+06  7.3547E-01 
.13800  2.025 1E+06  7.4430E-01 
.13900  2.0250E+06  7.5312E-01 
.14000  2.0250E+06  7.6194E-01 
.14100  2.0249E+06  7.7077E-01 
.14200  2.0248E+06  7.7959E-01 
.14300  2.0247E+06  7.8841E-01 
.14400  2.0246E+06  7.9724E-01 
.14500  2.0245E+06  8.0606E-01 
.14600  2.0244E+06  8.1488E-01 
.14700  2.0243E+06  8.2371E-01 
.14800  2.0242E+06  8.3253E-01 
.14900  2.024 1E+06  8.4135E-01 
.15000  2.0240E+06  8.5017E-01 
.15100  2.0239E+06  8.5900E-01 
.15200  2.023 8E+06  8.6782E-01 
.15300  2.0237E+06  8.7664E-01 
.15400  2.0236E+06  8.8546E-01 
.15500  2.023 5E+06  8.9428E-01 
.15600  2.0234E+06  9.0311E-01 
.15700  2.0233E+06  9.1193E-01 
.15800  2.0232E+06  9.2075E-01 
.15900  2.023 1E+06  9.2957E-01 
.16000  2.0230E+06  9.3839E-01 
.16100  2.0229E+06  9.4721E-01 
.16200  2.0227E+06  9.5603E-01 
.16300  2.0226E+06  9.6485E-01 
.16400  2.0225E+06  9.7367E-01 
.16500  2.0224E+06  9.8249E-01 
.16600  2.0223E+06  9.913  IE-01 


9.7835E-01 
9.7934E-01 
9.7982E-01 
9.8076E-01 
9.8168E-01 
9.8258E-01 
9.8345E-01 
9.8387E-01 
9.8471E-01 
9.8553E-01 
9.8632E-01 
9.8710E-01 
9.8748E-01 
9.8823E-01 
9.8896E-01 
9.8967E-01 
9.9037E-01 
9.9105E-01 
9.9138E-01 
9.9204E-01 
9.9268E-01 
9.933  IE-01 
9.9392E-01 
9.9423E-01 
9.9482E-01 
9.9540E-01 
9.9597E-01 
9.9653E-01 
9.9680E-01 
9.9734E-01 
9.9787E-01 
9.9839E-01 
9.9890E-01 
9.991 5E-01 
9.9964E-01 


2.1730E+01 

2.2008E+01 

2.2286E+01 

2.2564E+01 

2.2842E+01 

2.3121E+01 

2.3399E+01 

2.3677E+01 

2.3956E+01 

2.4234E+01 

2.4512E+01 

2.4791E+01 

2.5069E+01 

2.5348E+01 

2.5626E+01 

2.5905E+01 

2.6184E+01 

2.6462E+01 

2.6741E+01 

2.7020E+01 

2.7298E+01 

2.7577E+01 

2.7856E+01 

2.8135E+01 

2.8414E+01 

2.8692E+01 

2.8971E+01 

2.9250E+01 

2.9529E+01 

2.9808E+01 

3.0087E+01 

3.0366E+01 

3.0646E+01 

3.0925E+01 

3.1204E+01 


4.9780E+00 

5.0482E+00 

5.0833E+00 

5.1535E+00 

5.2237E+00 

5.2939E+00 

5.3641E+00 

5.3992E+00 

5.4694E+00 

5.5396E+00 

5.6098E+00 

5.6801E+00 

5.7152E+00 

5.7854E+00 

5.8556E+00 

5.9259E+00 

5.9961E+00 

6.0664E+00 

6.1015E+00 

6.1718E+00 

6.2420E+00 

6.3123E+00 

6.3826E+00 

6.4177E+00 

6.4880E+00 

6.5583E+00 

6.6286E+00 

6.6989E+00 

6.7341E+00 

6.8044E+00 

6.8747E+00 

6.9450E+00 

7.0153E+00 

7.0505E+00 

7.1208E+00 


Complete  Evaporation  at  z  =  .167 0(m) 


Superheated  Single  Phase  Results 

z  (m)  T  (K)  P  (Pa)  H  (J/kg)  (-) 

.16800  495.75  2.0220E+06  2.8147E+06 
.16900  505.06  2.0218E+06  2.8313E+06 
.17000  514.33  2.0216E+06  2.8480E+06 
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.17100  523.58  2.0214E+06  2.8647E+06 
.17200  530.76  2.0212E+06  2.8813E+06 
.17300  540.73  2.021 1E+06  2.8980E+06 
.17400  550.62  2.0209E+06  2.9146E+06 
.17500  567.63  2.0207E+06  2.9313E+06 
.17600  567.63  2.0205E+06  2.9479E+06 
.17700  573.42  2.0203E+06  2.9646E+06 
.17800  580.50  2.0201E+06  2.9813E+06 
.17900  588.76  2.0200E+06  2.9979E+06 
.18000  597.82  2.0198E+06  3.0146E+06 
.18100  610.65  2.0196E+06  3.0312E+06 
.18200  623.43  2.0194E+06  3.0479E+06 
.18300  623.43  2.0192E+06  3.0645E+06 
.18400  629.81  2.0190E+06  3.0812E+06 
.18500  639.49  2.0189E+06  3.0979E+06 
.18600  646.33  2.0187E+06  3.1145E+06 
.18700  659.97  2.0185E+06  3.1312E+06 
.18800  666.76  2.0183E+06  3.1478E+06 
.18900  673.54  2.0181E+06  3.1645E+06 
.19000  677.86  2.0179E+06  3.1812E+06 
.19100  686.94  2.0178E+06  3.1978E+06 
.19200  696.50  2.0176E+06  3.2145E+06 
.19300  703.30  2.0174E+06  3.2311E+06 
.19400  710.09  2.0172E+06  3.2478E+06 
.19500  723.63  2.0170E+06  3.2644E+06 
.19600  726.46  2.0169E+06  3.281 1E+06 
.19700  735.64  2.0167E+06  3.2978E+06 
.19800  748.35  2.0165E+06  3.3144E+06 
.19900  755.10  2.0163E+06  3.3311E+06 
.20000  764.57  2.0161E+06  3.3477E+06 
.20100  764.57  2.0161E+06  3.3477E+06 
.20200  764.57  2.0161E+06  3.3477E+06 
.20300  764.57  2.0161E+06  3.3477E+06 
.20400  764.57  2.0161E+06  3.3477E+06 
.20500  764.57  2.0160E+06  3.3477E+06 
.20600  764.57  2.0160E+06  3.3477E+06 
.20700  764.57  2.0160E+06  3.3477E+06 
.20800  764.57  2.0160E+06  3.3477E+06 
.20900  764.57  2.0160E+06  3.3477E+06 
.21000  764.57  2.0160E+06  3.3477E+06 


PEXIT=  .201E+07 
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Section  1 
Introduction 


1.1  Problem  Background 

DROPCON  is  a  computer  code  that  has  been  developed  to  model  the  condensation 
and  condensate  droplet  growth  phenomena  resulting  from  the  depressurization  of  a  vapor- 
noncondensable  mixture.  The  code  is  meant  to  be  a  generic  analytical  tool  for  the 
calculation  of  exit  conditions  of  components  such  as  turbine  expanders. 

The  code  assumes  that  a  known  flow  rate  of  a  vapor-noncondensable  mixture,  initially  at 
saturation  conditions,  flows  into  a  device  that  imposes  a  pressure  drop  to  the  flow  field. 
The  conditions  of  the  flow  mixture  at  the  exit  of  the  device  are  then  calculated  by  the 
code.  Parameters  that  are  calculated  are  meant  to  be  sufficient  for  modeling  the  droplet 
removal  phenomena  in  a  separator/dryer  component. 

1.2  Report  Summary 

Section  2  briefly  describes  the  phenomena  that  are  modeled  in  this  computer  code.  The 
code  structure  along  with  its  flow  diagram  is  presented  in  Section  3.  Section  4  briefly 
describes  each  of  the  routines  used  in  the  code.  The  input  and  output  are  described  in 
Section  5.  A  sample  case  is  provided  in  Section  6,  where  a  complete  listing  of  all  the 
input  and  output  is  included.  The  references  are  listed  in  Section  7. 
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Section  2 

Description  of  Theory 


A  schematic  diagram  of  the  modeled  system  is  shown  in  Figure  1.  The  system  is  basically 
a  control  volume  with  an  arbitrary  configuration.  Furthermore,  it  is  assumed  to  exchange 
mechanical  work  with  its  surroundings.  At  Port  1,  the  inlet  port,  a  known  flow  rate  of  a 
saturated  mixture  of  a  vapor  (steam  at  this  time)  and  a  non-condensable  (air  at  this  time) 
with  known  thermodynamic  conditions  (pressure,  vapor  mass  fraction)  is  mixed  with  a 
known  flow  rate  of  saturated  droplets  sprayed  into  the  flow  field.  The  mixture  is  then 
assumed  to  undergo  depressurization,  and  leave  the  systems  at  Port  2,  the  exit  port,  while 
it  has  reached  at  thermodynamic  equilibrium.  It  is  also  assumed  that  the  pressure  drop  in 
the  system  is  known. 

The  model  is  based  on  the  following  assumptions. 

1 .  Thermodynamic  equilibrium  is  maintained  throughout  the  system. 

2.  Liquid  is  impermeable  to  the  non-condensable. 

3.  The  non-condensable  is  an  ideal  gas. 

4.  Excluding  the  inlet  and  outlet  ports,  no  mass  transfer  takes  place  anywhere  else 
between  the  modeled  system  and  the  surroundings. 

5.  The  liquid-gas  mixture  is  homogeneous  at  both  inlet  and  exit  ports. 

6.  If  droplets  are  injected  into  the  stream  at  inlet,  the  injected  spray  droplets  at  inlet  have 
a  log-normal  size  distribution  with  known  distribution  characteristics. 

Consistent  with  the  above  assumptions,  the  conservation  equations  for  the  total  mass,  the 
condensable  (liquid  +  vapor)  mass,  and  total  energy  (thermal  and  mechanical)  are 
iteratively  solved,  thereby  the  exit  thermodynamic  conditions,  as  well  as  the  mass  flow 
rates  of  vapor,  non-condensable  and  droplets,  are  calculated. 

Four  different  types  of  expanders  are  modeled.  The  following  table  is  a  summary  of  the 
properties  of  the  modeled  expander  types. 
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Expander  Type 

Description 

1 

Inlet  and  outlet  pressures  are  known;  droplets  with  known  size 
distribution  are  injected  at  inlet. 

2 

Turbine  with  known  Input  pressure,  output  power,  and  isentropic 
efficiency;  fog  formation  takes  place. 

3 

Turbine  with  known  Input  pressure,  output  power,  and  isentropic 
efficiency;  droplets  with  known  size  distribution  are  injected  at 
inlet. 

4 

Inlet  and  outlet  pressures  are  known;  fog  formation  takes  place. 

The  injected  droplets,  when  applicable,  are  assumed  log-normally  distributed  with  respect 
to  their  diameters.  The  droplets  leaving  the  system  are  also  assumed  to  have  a  log¬ 
normal  size  distribution,  with  a  geometric  standard  deviation  equal  to  that  at  the  inlet. 

This  is  equivalent  to  the  assumption  that  the  shape  of  the  size  distribution  does  not 
change,  and  only  the  distribution  curve  moves  towards  a  lager  size.  The  code  thus 
calculates  the  size  distribution  characteristics  of  the  droplets  at  the  exit. 

At  present  the  code  is  equipped  with  property  routine  for  water/steam  (the  condensate) 
and  air  (the  non-condensable). 
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Figure  1 

Schematic  of  the  Modeled  System 


Section  4 


Description  of  Subroutines  and  Functions 

This  section  contains  a  description  of  subroutines  and  functions  used  in  this  computer 
code. 

Subroutines: 

TURVl  Models  the  first  type  of  turbine  expander. 

FS1  Calculates  the  function  values  for  the  first  type  turbine. 

TURV2  Models  the  second  type  of  turbine  expander. 

FS2  Calculates  the  function  values  for  the  second  type  turbine. 

TURV3  Models  the  third  type  of  turbine  expander. 

COE  Calculates  the  coefficients  of  each  equation  at  each  iteration  for 

the  FS1,FS2  and  FS3  subroutines. 

FS3  Calculates  the  function  values  for  the  third  type  turbine. 

TURV4  Models  the  fourth  type  of  turbine  expander. 

FS4  Calculates  the  function  values  for  the  fourth  type  turbine. 

PNEWTON  Solves  the  set  of  non-linear  equations. 

AGAUS  Solves  the  set  of  linear  equations. 

POLINT  This  is  a  Lagrangian  interpolation  routine,  is  used  to  interpolate 
between  thermodynamic  values  in  property  tables. 
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Functions: 

TSATF 

Given  the  pressure,  this  returns  the  saturation  temperature.  T  in  K, 
P  in  bar. 

VFT 

Given  the  saturation  temperature,  this  finds  the  density  of  water.  T 
in  K,  density  in  kg/m3. 

HFPSAT 

Finds  the  liquid  saturation  enthalpy  given  the  pressure.  P  in  bar, 
enthalpy  in  J/kg. 

HGPSAT 

Finds  the  vapor  saturation  enthalpy  given  the  pressure.  P  in  bar, 
enthalpy  in  J/kg. 

HAT29 

Finds  enthalpy  for  non-condensible  AIR  given  the  temperature.  T 
in  K,  enthalpy  in  kJ/kg. 

VGPT 

Finds  specific  volume  given  pressure  and  temperature,  P  in  bar,  T 
in  K,  and  VG  in  cm3/g. 

SAT29 

Finds  entropy  of  non-condensible  AIR  given  the  temperature.  T 
in  K  and  SAT29  in  J/kg.K. 

SGPSAT 

Finds  entropy  of  saturated  vapor  in  partial  pressure  of  vapor  given 
the  temperature.  T  in  C,  and  entropy  in  J/kg.K 

SFPSAT 

Finds  entropy  of  saturated  liquid  (water)  in  vapor  partial  pressure 
given  the  temperature.  T  in  C,  and  entropy  in  J/kg.K 
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Section  5 


Input  and  output 


This  section  describes  the  input  parameters  necessary  to  run  this  computer  code  and  an 
explanation  of  all  the  associated  output. 


INPUT: 

Table  2  defines  all  the  input  parameters  for  TURV1,  the  first  expander  type.  For  other 
expander  types  the  definition  of  parameters  are  the  same.  The  sample  input  data  files,  to 
be  given  shortly,  provide  the  exact  order  of  input  parameter  that  should  be  provided. 


Table  2 

Parameters  units  Description 


w 

J 

Turbine  power 

PI 

bar 

Inlet  pressure 

P2 

bar 

Outlet  pressure 

MF1 

kg/s 

Mass  flow  rate  of  liquid  at  inlet 

El 

Volume  fraction  occupied  by  liquid 

NT1 

m'3 

Total  number  density  of  droplets  at  inlet 

SIG 

Geometric  standard  deviation 

R 

m3.bar/(kmol.K) 

Universal  gas  constant 

MN 

kg/kmole 

Molar  mass  of  non-condensable 

MV 

kg/kmole 

Molar  mass  of  water 

U1 

m/s 

Mixture  velocity  at  inlet 

A2 

m2 

Outlet  area 

MG1 

kg/s 

Mass  flow  rate  of  gas  mixture  at  inlet 

MV1 

kg/s 

Mass  flow  rate  of  vapor  at  inlet 

MN1 

kg/s 

Mass  flow  rate  of  non-condensable  at  inlet 

PREF 

bar 

Reference  pressure  for  entropy 

TREF 

K 

Reference  temperature  for  entropy 
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OUTPUT: 


Two  different  types  of  output  are  generated  by  the  code.  One  type  of  output  file,  with  the 
generic  name  TURVn.out,  where  n  represents  the  expander  type,  provides  detailed  and 
well-defined  information  regarding  the  exit  conditions  of  the  expander.  Another  output 
file  with  generic  name  CYC-n.dat,  with  n  representing  the  expander  type,  is  also 
generated.  This  file  can  be  used  as  the  input  file  for  the  code  CYCSEP  which  models  a 
cyclone-type  separator.  Only  the  parameters  in  the  former  data  files  are  described  below. 

Table  3  defines  all  output  parameters  for  TURVl.out,  TURV2.out,  TURV3.out,  and 
TURV4.out.  (They  are  all  formatted  the  same  way.) 


Table  3 


Parameters 

Units 

Description 

U2 

m/s 

Outlet  velocity  of  gas  mixture 

A2 

m 

Outlet  area 

W 

j 

turbine  power 

E2 

Volume  fraction  occupied  by  liquid 

XN2 

Mass  fraction  of  noncondensable  in  gas  mixture 

YN2 

Mole  fraction  of  noncondensable  in  gas  mixture 

PV2 

bar 

Partial  pressure  of  vapor  in  gas  mixture 

PN2 

bar 

Partial  pressure  of  noncondensable  in  gas  mixture 

T2 

K 

Temperature  of  gas  mixture 

HG2 

J/kg 

Enthalpy  of  gas  mixture 

RHOG2 

kg/m3 

Density  of  gas  mixture 

SIG 

Geometric  standard  deviation 

INCMD 

Log  median  diameter 

DV 

m 

Volume-averaged  droplet  diameter 

10 


Section  6 

Sample  Input  and  Output 


6.1.1  Sample  Input  For  the  First  Expander  Type 


0.0,5. 5,  1.5,0.02, 

l.e-6, l.e08, 2. 0,0. 0831434, 29. 0,18. 0,2. 0,2. 02683E-3, 
0.08, 0.04, 0.04, 

1.0,100.0 


6.1.2  Sample  Output  For  the  First  Expander  Type 


OUTLET: 

U2=  .362170D+02 
A 2=  .202683D-02 
W=  .000000D+00 
E2=  .286907D-03 
XN2=  .501698D+00 
YN2=  .384585D+00 
PV2=  .923123D+00 
PN2=  .576877D+00 
T2=  .370557D+03 
HG2=  .151829D+07 
RHOG2=  .108646D+01 

GEOMETRIC  STANDARD  DEVIATION  IS  .200000D+01 
LOG  MEDIAN  DIAMETER  IS:  -.100847D+02 
MEDIAN  OF  VOLUME  DISTRIBUTION  IS:  .176298D-03 
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6.2.1  Sample  Input  For  the  Second  Expander  Type 


5696.46,5.5,  1.3, 

l.e-6,l.e08, 2.0, 0.0831434, 29.0, 18.0, 2.0, 2.02683E-3, 
0.08,0.04,0.04, 

1.0,100.0 


6.2.2  Sample  Input  For  the  Second  Expander  Type 

OUTLET: 

U2=  .405622D+02 
A2=  .202683D-02 
W=  .569646D+04 
E2=  .184247D-04 
XN2=  .509310D+00 
YN2=  .391817D+00 
PV2=  .790637D+00 
PN2=  .509363D+00 
T2=  .366344D+03 
HG2=  .149590D+07 
RHOG2=  .955316D+00 
SG2=  .349689D+04 
SF2=  .219397D+04 

GEOMETRIC  STANDARD  DEVIATION  IS  .200000D+0 1 
LOG  MEDIAN  DIAMETER  IS:  -.1 09998D+02 
MEDIAN  OF  VOLUME  DISTRIBUTION  IS:  .705994D-04 
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6.3.1  Sample  Input  For  the  Third  Expander  Type 


3580.12,5.5,  1.3,0.02, 

1  .e-6, 1  ,e08, 2.0, 0.083 1434,29.0, 1 8.0,2.0,2.02683E-3, 
0.08,0.04,0.04, 

1.0,100.0 


6.3.2  Sample  Input  For  the  Third  Expander  Type 

OUTLET: 

U2=  .406802D+02 
A2=  .202683D-02 
W=  .358012D+04 
E2=  .161475D-04 
XN2=  .508165D+00 
YN2=  .390726D+00 
PV2=  .792057D+00 
PN2=  .507943D+00 
T2=  .366392D+03 
HG2=  .149859D+07 
RHOG2=  .954689D+00 
SG2=  .349666D+04 
SF2=  .219397D+04 

GEOMETRIC  STANDARD  DEVIATION  IS  .200000D+01 
LOG  MEDIAN  DIAMETER  IS :  -.11043  8D+02 
MEDIAN  OF  VOLUME  DISTRIBUTION  IS:  .675620D-04 
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6.4.1  Sample  Input  For  the  Fourth  Expander  Type 


0.0,5 .5,  2.5, 

1  .e-6, 1  ,e08, 2.0, 0.083 1434,29.0,1 8.0, 2.0, 2.02683E-3, 
0.08,0.04,0.04, 

1.0,100.0 


6.4.2  Sample  Input  For  the  Fourth  Expander  Type 

OUTLET: 

U2=  .225896D+02 
A2=  .202683D-02 
W=  .OOOOOOD+OO 
E2=  .299214D-05 
XN2=  .500816D+00 
YN2=  .383750D+00 
PV2=  .154062D+01 
PN2=  .959376D+00 
T2=  .385332D+03 
HG2=  .152803D+07 
RHOG2=  .174445D+01 

GEOMETRIC  STANDARD  DEVIATION  IS  .200000D+01 
LOG  MEDIAN  DIAMETER  IS :  -.11 6057D+02 
MEDIAN  OF  VOLUME  DISTRIBUTION  IS:  .385178D-04 
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Section  1 

1.  Introduction 


1.1  Problem  Background 

CYCSEP  is  a  computer  code  that  has  been  developed  to  model  the  condensation  and 
droplet  removal  in  a  system  consisting  of  a  cyclone-type  separator.  Condensation  takes 
place  as  a  result  of  the  depressurization  of  a  vapor-noncondensable  mixture,  while  droplet 
removal  is  caused  by  the  impaction  of  droplets  with  the  cyclone  walls. 

CYCSEP  should  be  used  in  tandem  with  the  computer  code  DROPCON,  and  uses  two 
input  files.  One  input  file,  named  CYCGMTRY.DAT,  contains  geometrical  parameters 
for  the  cyclone  separator.  In  addition,  CYCSEP  uses  an  output  file  directly  generated  by 
the  DRPCOND  code  package.  Therefore,  before  CYCSEP  is  used,  DROPCON  must  be 
run,  so  that  the  aforementioned  input  data  file  for  CYCSEP  is  automatically  generated. 
Thus,  a  review  of  the  user’s  manual  for  DROPCON  is  necessary  before  the  present  code  is 
used. 

Figure  1  is  a  schematic  of  the  modeled  separator.  The  separator  model  is  a 
modification  of  a  classical  and  widely-applied  model  (Leith  and  Licht,  1972).  The  latter 
model,  in  its  original  form,  applies  to  the  Figure  1(a).  The  modified  model,  however,  is 
meant  to  represent  an  experimental  test  section  displayed  in  Figure  1(b),  and  currently 
under  development  at  Georgia  Tech. 

1.2  Report  Summary 

The  phenomena  that  are  modeled  in  this  computer  code  package  are  described  in 
Ghiaasiaan  et  al  (2002)  and  will  not  be  repeated  here.  Section  2  briefly  describes  each  of 
the  routines  used  in  the  code  package.  The  input  and  output  are  described  in  Section  3. 

A  sample  case  is  provided  in  Section  4,  where  a  complete  listing  of  all  the  output  is 
included.  The  references  are  listed  in  Section  5. 
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Section  2 

Description  of  Subroutines  and  Functions 

This  section  contains  a  description  of  each  subroutine  and  function  in  this  computer  code. 
The  reader  is  reminded  that  this  code  is  meant  to  use  the  output  provided  by  the  code 
package  DROPCON.  The  latter  code  package  contains  models  for  four  different 
expanders.  The  expander  types  referred  to  in  the  forthcoming  list  are  briefly  described  in 
the  DROPCON  manual,  and  in  Ghiaasiaan  et  al.  (2002).  These  documents  should  be 
consulted  before  the  present  code  is  applied. 


Subroutines: 


CYC1 

CYC2 

CYC3 

CYC4 

COE 

FS 


Calculates  the  cyclone  results  when  the  cyclone  is  attached  to  the 
first  turbine  type. 

Calculates  the  cyclone  results  when  the  cyclone  is  attached  to  the 
second  turbine  type. 

Calculates  the  cyclone  results  when  the  cyclone  is  attached  to  the 
third  turbine  type. 

Calculates  the  cyclone  results  when  the  cyclone  is  attached  to  the 
fourth  turbine  type. 

Calculates  the  coefficients  of  each  equation  to  be  solved  at  each 
iteration. 

Calculates  the  numerical  value  of  each  equation  using  the  results  of 
COE. 


PNEWTON  Solves  the  set  of  non-linear  equations  using  Newton’s  method. 


AGAUS 

POLINT 

VISCOV 


Solves  the  set  of  linear  equations. 

This  is  a  Lagrangian  interpolation  routine,  and  is  used  to  interpolate 
between  thermodynamic  values  in  property  tables. 

Gets  vapor  viscosity  EMU  at  given  pressure  P  and  temperature 
T,  P  in  bar,  T  in  K,  EMU  in  Pa.s. 


Functions: 

TSATF  Given  the  pressure,  this  returns  the  saturation  temperature.  T  in  K, 

P  in  bar. 
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VFT 

Given  the  saturation  temperature,  this  finds  the  density  of  water.  T 
in  K,  density  in  kg/m3. 

HFPSAT 

Finds  the  liquid  saturation  enthalpy  given  the  pressure.  P  in  bar, 
enthalpy  in  J/kg. 

HGPSAT 

Finds  the  vapor  saturation  enthalpy  given  the  pressure.  P  in  bar, 
enthalpy  in  J/kg. 

HAT29 

Finds  enthalpy  for  non-condensible  AIR  given  the  temperature.  T I 
In  K,  enthalpy  in  kJ/kg. 

VGPT 

Finds  specific  volume  given  pressure  and  temperature,  P  in  bar,  T 
in  K,  and  VG  in  cm3/g. 
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Section  3 
Input  and  output 


This  section  describes  the  input  parameters  necessary  to  run  this  computer  code,  and 
provides  an  explanation  of  all  the  associated  output. 

INPUT 

The  input  files  for  this  code,  as  mentioned  earlier,  include  two  separate  files.  One  input 
file,  named  CYCGMTRY.DAT,  contains  geometric  parameters  and  outlet  pressure.  In 
addition  another  input  file  is  automatically  generated  by  the  code  package  DRPCOND. 
Since  the  code  uses  input  provided  by  four  different  expanders,  the  user  must  first  specify 
the  expander  type.  (For  expander  types,  see  the  User’s  Manual  for  DROPCON  Code). 

Two  different  types  of  output  files  are  generated  by  the  DROPCON  code.  One  type  of 
output  file,  with  the  generic  name  TURVn.out,  where  n  represents  the  expander  type, 
provides  detailed  and  well-defined  information  regarding  the  exit  conditions  of  the 
expander.  Another  output  file  with  generic  name  CYC-n.dat,  with  n  representing  the 
expander  type,  is  also  generated  by  DROPCON.  This  file  is  meant  to  be  used  as  one  of 
the  input  files  for  the  code  CYCSEP  code. 

The  specification  of  the  expander  type  by  the  user  is  done  interactively.  The  CYCSEP 
code,  once  executed,  enquires  on  the  screen: 

“Please  enter  the  expander  type  (1-4)?” 

In  response,  the  user  must  type  the  expander  type  number,  followed  by  a  return.  The 
code  then  does  the  following.  If  expander  type  1  is  specified  by  the  user,  the  code 
searches  for  CYC-l.dat,  and  reads  it.  In  case  expander  type  is  2,  then  the  code  searches 
for  CYC-2.dat,  and  so  on. 

Transfer  of  information  from  DROPCON  to  CYCSEP  is  thus  automatic,  as  long  as  the 
correct  expander  type  is  specified  by  the  user.  A  list  of  parameters  provided  by  CYC- 
n.dat  files  will  therefore  not  be  provided  here. 

Table  1  is  the  list  of  parameters  that  are  provided  in  the  CYCGMTRY.DAT  data  file. 
Numerical  values  of  these  parameters  must  be  typed  into  the  file  CYCGMTRY.DAT,  in 
the  same  order  as  they  appear  in  the  table,  with  numbers  separated  from  each  other  by  a 
comma. 

Table  2  is  a  list  of  the  parameters  that  appear  in  the  output  from  CYCSEP  execution. 
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Table  1 

List  of  Parameters  Provided  in  CYCGMTRY.DAT 


Parameters  units  Description 


D 

m 

Body  diameter  of  cyclone 

DE 

m 

Pipe  diameter  at  inlet  and  outlet 

DF 

ft 

Diameter  of  cyclone  body 

HI 

m 

Height  of  cylindrical  body 

H2 

m 

Total  cyclone  height 

VS 

rv.3 

m 

Annular-shaped  volume  above  exit  duct  to  mid  level  of 
entrance  duct 

MN 

kg/kmole 

Molar  mass  of  noncondensable 

MV 

kg/kmole 

Molar  mass  of  water  vapor 

P2 

bar 

Mixture  pressure  at  outlet 

SIG 

Geometric  standard  deviation 
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Table  2 

List  of  Parameters  in  the  Output  Files 


Parameters 

Units 

Description 

A2 

m 

Area  of  outlet 

U2 

m/s 

Velocity  of  mixture 

E2 

Volume  fraction  occupied  by  liquid 

XN2 

Mass  fraction  of  noncondensable  in  gas  mixture 

YN2 

Mole  fraction  of  noncondensable  in  gas  mixture 

PV2 

bar 

Partial  pressure  of  vapor  in  gas  mixture 

PN2 

bar 

Partial  pressure  of  noncondensable  in  gas  mixture 

T2 

K 

Temperature  of  gas  mixture 

RHOF2 

kg/m3 

Density  of  droplets 

RHOG2 

kg/m3 

Density  of  gas  mixture 

HG2 

J/kg 

Enthalpy  of  gas  mixture 

HF2 

J/kg 

Enthalpy  of  droplets 

EMUG2 

Pa.S 

Viscosity  of  gas  mixture 

SIG 

Geometric  standard  deviation 

INCMD 

Log  mean  diameter 

DV2 

m 

droplets  volumetric  meandiameter 

NT2 

#/m3 

Total  number  density  of  droplets  at  outlet 

ETA 

Droplet  removal  efficiency 
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Section  4 

Sample  Input  and  Output 


Sample  input  and  output  files  for  all  options  are  provided  in  this  section.  Note  that  the 
input  file  CYCGMTRY.DAT  is  common  for  all  of  the  options.  Also,  note  that  files  with 
generic  name  cyc-n.dat  are  generated  by  the  DROPCON  computer  code. 

4.1.1  Sample  CYCGMTRY.DAT  Input  File 

0.168275, 5.08E-002, 0.1008,0.4318,0.53848, 1.3E-003 
29.0,18.0,1.0,2.0 

4.1.2  Sample  Input  cyc-l.dat 

0.000000000000000E+000  2.026829769984279E-003  36.216971074700500 

2 . 8690726600263 86E-004  5.016977999359462E-001  3.84584801 1661 990E-001 
9.231227982507015E-001  5.768772017492985E-001  370.557208463130800 

1 . 000000000000000E+008  1 .76298 1 5 1 3684647E-004 
1518290.426654123000000 
1.086458155699219 

4.1.3  Sample  Output  cyc-l.out 

THE  FOLLOWING  IS  THE  CYCLONE  RESULT  FROM  THE  TURBINE 
FIRST  TYPE,  W=0  MF1>0,  OUTPUT  DATA 
OUTLET: 

A2=  .202683D-02 
U2=  .535646D+02 
E2=  .37501 8D-1 8 
XN2=  .494655D+00 
YN2=  .377939D+00 
PV2=  . 62206 1D+00 
PN2=  .377939D+00 
T2=  .360030D+03 
RHOF2=  .969690D+03 
RHOG2=  .742217D+00 
HG2=  .152085D+07 
HF2=  .368708D+06 
EMUG2=  .152999D-04 

GEOMETRIC  STANDARD  DEVIATION  IS  .200000D+0 1 
LOG  MEDIAN  DIAMETER  IS:  -.10101 1D+02 
MEDIAN  OF  VOLUME  DISTRIBUTION  IS :  .1 73435D-03 

TOTAL  DROPLET  NUMBER  AT  OUTLET  IS:  .  144329D-06 
REMOVAL  EFFICIENCY  ETA  IS:  .  100000D+01 
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4.2.1  Sample  Input  cyc-2.dat 


5696.456450794689000  2.026829769984279E-003  40.562175665823850 

1.842470350576125E-005  5.093 101708643396E-001  3.918174398295731E-001 
7.90637299221 1506E-001  5.093626530951335E-001  366.343828883619500 

1.000000000000000E+008  7.059935695044410E-005  1495899.015409790000000 
9.553 161467587280E-001 


4.2.2  Sample  Output  cyc-2.out 

THE  FOLLOWING  IS  THE  CYCLONE  RESULT  FROM  THE  TURBINE 
SECOND  TYPE,  W>0  MFL=0,  OUTPUT  DATA 
OUTLET: 

A2=  .202683D-02 
U2=  .520967D+02 
E2=  .118021D-12 
XN2=  .505579D+00 
YN2=  .388266D+00 
PV2=  .611734D+00 
PN2=  .388266D+00 
T2=  .359598D+03 
RHOF2=  .969973D+03 
RHOG2=  .746845D+00 
HG2=  .149514D+07 
HF2=  .366588D+06 
EMUG2=  .153796D-04 

GEOMETRIC  STANDARD  DEVIATION  IS  .200000D+01 
LOG  MEDIAN  DIAMETER  IS:  -.1 10841D+02 
MEDIAN  OF  VOLUME  DISTRIBUTION  IS:  .648970D-04 
TOTAL  DROPLET  NUMBER  AT  OUTLET  IS:  .108665D+01 
REMOVAL  EFFICIENCY  ETA  IS:  .100000D+01 


4.3.1  Sample  Input  cyc-3.dat 


5114.459519643379000  2.026829769984279E-003  40.680241677380760 

1.614753308588079E-005  5.081645586989505E-001  3.907256674766556E-001 
7.920566032278834E-001  5.079433490884006E-001  366.391954182002000 

1.000000000000000E+008  6.7562021 74997793E-005  1498594.811873642000000 
9.546887900095444E-001 


4.3.2  Sample  Output  cyc-3.out 

THE  FOLLOWING  IS  THE  CYCLONE  RESULT  FROM  THE  TURBINE 
THIRD  TYPE,W>0  MFL>0,OUTP  UT  DATA 
OUTLET: 

A2=  .202683D-02 
U2=  .522460D+02 
E2=  .184953D-12 
XN2=  .50445 1D+00 
YN2=  .387195D+00 
PV2=  .612805D+00 
PN2=  .387195D+00 
T2=  .359643D+03 
RHOF2=  .969944D+03 
RHOG2=  .746364D+00 
HG2=  .149779D+07 
HF2=  .366810D+06 
EMUG2=  .153713D-04 

GEOMETRIC  STANDARD  DEVIATION  IS  .200000D+01 
LOG  MEDIAN  DIAMETER  IS :  -.111 406D+02 
MEDIAN  OF  VOLUME  DISTRIBUTION  IS:  .613306D-04 
TOTAL  DROPLET  NUMBER  AT  OUTLET  IS:  .210937D+01 
REMOVAL  EFFICIENCY  ETA  IS:  .100000D+01 
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4.4.1  Sample  Input  cyc-4.dat 


O.OOOOOOOOOOOOOOOE+OOO  2.026829769984279E-003  22.589572022357500 

2.992 1 43576350699E-006  5 .008 1 60082856482E-001  3 .837503288 1 14426E-001 
1.540624177971394  9.593758220286065E-001  385.331703469266300 

1.000000000000000E+008  3.851779218993663E-005  1528026.628694685000000 
1.744447088390374 


4.4.2  Sample  Output  cyc-4.out 

THE  FOLLOWING  IS  THE  CYCLONE  RESULT  FROM  THE 

TURBINE  FOURTH  TYPE,  W=0  MFL=0, 0 

UTPUT  DATA 

OUTLET: 

A2=  .202683D-02 
U2=  .535646D+02 
E2=  .000000D+00 
XN2=  .494655D+00 
YN2=  .377939D+00 
PV2—  .622061D+00 
PN2=  .377939D+00 
T2=  .360030D+03 
RHOF2=  .969690D+03 
RHOG2=  .742217D+00 
HG2=  .152085D+07 
HF2=  .368708D+06 
EMUG2=  .152999D-04 

GEOMETRIC  STANDARD  DEVIATION  IS  .200000D+01 
LOG  MEDIAN  DIAMETER  IS:  -.  107778D+02 
MEDIAN  OF  VOLUME  DISTRIBUTION  IS:  .881491D-04 
TOTAL  DROPLET  NUMBER  AT  OUTLET  IS:  .144329D-06 
REMOVAL  EFFICIENCY  ETA  IS:  .100000D+01 
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